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Preface 



These notes are the written version of lectures given at the Compstar 2010 School held in Caen 
on 8-16 February 2010. They are intended to introduce the theory of rotating stars in general 
relativity. Whereas the application is clearly towards neutron stars (and associated objects 
like strange quark stars), these last ones are not discussed here, except for illustrative purposes. 
Instead the focus is put on the theoretical foundations, with a detailed discussion of the spacetime 
symmetries, the choice of coordinates and the derivation of the equations of structure from the 
Einstein equation. The global properties of rotating stars (mass, angular momentum, redshifts, 
orbits, etc.) are also introduced. These notes can be fruitfully complemented by N. Stergioulas' 
review [95] and by J.L. Friedman & N. Stergioulas' textbook [44]. In addition, the textbook 
by P. Haensel, A. Y. Potekhin &; D. G. Yakovlev [55] provides a good introduction to rotating 
neutron star models computed by the method exposed here. 

The present notes are limited to perfect fluid stars. For more complicated physics (super- 
fluidity, dissipation, etc.), see N. Andersson & G. L. Comer's review [5]. Magnetized stars are 
treated in J. A. Pons' lectures [87]. Moreover, the stability of rotating stars is not investigated 
here. For a review of this rich topic, see N. Stergioulas' article [95] or L. Villain's one [103]. 

To make the exposure rather self-contained, a brief overview of general relativity and the 
3+1 formalism is given in Chap. 1. For the experienced reader, this can be considered as a mere 
presentation of the notations used in the text. Chap. 2 is entirely devoted to the spacetime 
symmetries relevant for rotating stars, namely stationarity and axisymmetry. The choice of co- 
ordinates is also discussed in this chapter. The system of partial differential equations governing 
the structure of rotating stars is derived from the Einstein equation in Chap. 3, where its numer- 
ical resolution is also discussed. Finally, Chap. 4 presents the global properties of rotating stars, 
some of them being directly measurable by a distant observer, like the redshifts. The circular 
orbits around the star are also discussed in this chapter. 
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1.1 Geometrical framework 

1.1.1 The spacetime of general relativity 

Relativity has performed the fusion of space and iime, two notions which were completely distinct 
in Newtonian mechanics. This gave rise to the concept of spacetime^ on which both the special 
and general theory of relativity are based. Although this is not particularly fruitful (except for 
contrasting with the relativistic case), one may also speak of spacetime in the Newtonian frame- 
work: the Newtonian spacetime ^ is nothing but the affine space M^, foliated by the hyperplanes 
Tit of constant absolute time t: these hyperplanes represent the ordinary 3-dimensional space at 
successive instants. The foliation (Sf)tgig is a basic structure of the Newtonian spacetime and 
does not depend upon any observer. The worldline C oi a particle is the curve in ^ generated 
by the successive positions of the particle. At any point A E: the time read on a clock moving 
along C is simply the parameter t of the hyperplane Sf that intersects C at A. 

The spacetime ^ of special relativity is the same mathematical space as the Newtonian one, 
i.e. the affine space M^. The major difference with the Newtonian case is that there does not 
exist any privileged foliation {Tjt)t&M.- Physically this means that the notion of absolute time 
is absent in special relativity. However ^ is still endowed with some absolute structure: the 
metric tensor g and the associated light cones. The metric tensor is a symmetric bilinear 
form g on which defines the scalar product of vectors. The null (isotropic) directions of g 
give the worldlines of photons (the light cones). Therefore these worldlines depend only on the 
absolute structure g and not, for instance, on the observer who emits the photon. 

The spacetime ^ of general relativity differs from both Newtonian and special relativistic 
spacetimes, in so far as it is no longer the affine space but a more general mathematical 
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Figure 1.1: The vectors at two points p and q on the spacetime manifold ^ belong to two different 
vector spaces: the tangent spaces and 

structure, namely a manifold. A manifold of dimension 4 is a topological space such that 
around each point there exists a neighbourhood which is homeomorphic to an open subset of M'^. 
This simply means that, locally, one can label the points of ^ in a continuous way by 4 real 
numbers (a^")oe{o,i,2,3} (which are called coordinates). To cover the full several different 
coordinate patches {charts in mathematical jargon) can be required. 

Within the manifold structure the definition of vectors is not as trivial as within the affine 
structure of the Newtonian and special relativistic spacetimes. Indeed, only infinitesimal vectors 
connecting two infinitely close points can be defined a priori on a manifold. At a given point 
p G the set of such vectors generates a 4-dimensional vector space, which is called the 
tangent space at the point p and is denoted by £^p{^). The situation is therefore different 
from the Newtonian or special relativistic one, for which the very definition of an affine space 
provides a unique global vector space. On a manifold there are as many vector spaces as points 
p (cf. Fig. 1.1). 

Given a coordinate system (x") around some point p G one can define a vector basis 
(9a)aG{0,i,2,3} of the tangent space S^p{^) by asking that the components of the infinitesimal 
vector p^' that connect neighbouring points p (coordinates (x")) and p' (coordinates (x" +dx°')) 
are the infinitesimal increments dx" in the coordinates when moving from p to p': 

3 

p^' = "^dx"da = dx''da. (1.1) 
a=0 

For the second equality in the above equation, we have used Einstein summation convention 
on repeated indices, which amounts to omitting the S sign. We shall make a systematic use of 
this convention in these notes. The vector basis (da) is called the natural basis associated 
with the coordinates (x"). The notation da is inspired by d/dx"' and is a reminiscence of the 
intrinsic definition of vectors on a manifold as differential operators acting on scalar fields (see 
e.g. [104]). 

For any vector v G =y^(^) (not necessarily an infinitesimal one), one defines the components 
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(f") of V in the basis (da) according to the expansion 

v = v"da. (1.2) 

1.1.2 Linear forms 

Besides vectors, important objects on the manifold ^ are hnear forms. At each point p € 
a linear form is a an apphcation^ 

that is hnear: (cj, Xv + u) = X{uJ, v) + {uj, u). The set of aU hnear forms at p constitutes a 4- 
dimensional vector space, which is cahed the dual space of J^p{^) and denoted by ^p{^)* . A 
field of linear forms on ^ is called a 1-form. Given the natural basis {da) of S^p{^) associated 
with the coordinates there is a unique basis of denoted by (dx°), such that 



/3 



^1.4) 



where 5"^ is the Kronecker symbol : 5"^ = 1 if a = /? and otherwise. The basis (dx") 

is called the dual basis of the basis (da)- The notation (dx°) stems from the fact that if we 
apply the linear form dx" to the infinitesimal displacement vector (1.1), we get nothing but the 
number dx": 

(dx°,^') = (dx", dx^ dp) = dxl^idx", dp) = dx^ 6"^ = dx". (1.5) 
The dual basis can be used to expand any linear form cj, thereby defining its components oja- 

u=ujadx°'. (1.6) 

In terms of components, the action of a linear form on a vector takes then a very simple form: 



{u,v) =u;av'' . (1.7) 



Given a smooth scalar field on namely a smooth application / : ^ — t- M, there exists a 
1-form canonically associated with it, called the gradient of f and denoted df. It is the unique 
linear form which, once applied to the infinitesimal displacement vector pf)', gives the change in 
/ between points p and p': 

df := f{p') - f{p) = dfiW). (1.8) 

Since df = df /dx'^ dx", formula (1.5) implies that the components of the gradient with respect 
to the dual basis are nothing but the partial derivatives of / with respect to the coordinates (x°) 

d/ = ^dx". (1.9) 

5x" ^ ' 

Remark : In non-relativistic physics, the gradient is very often considered as a vector and not as a 
1-form. This is so because one associates implicitly a vector u to any 1-form lo via the Euclidean 
scalar product of R'^: Vw G R'^, (w, v) = lj-v. Accordingly, formula (1.8) is rewritten as df = V f-pf) . 
But one shall keep in mind that, fundamentally, the gradient is a linear form and not a vector. 



^We are using the same bra-ket notation as in quantum mechanics to denote the action of a linear form on a 
vector. 
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1.1.3 Tensors 

In relativistic physics, an abundant use is made of linear forms and their generahzations: the 
tensors. A tensor of type {k,£), also called tensor k times contravariant and i times 
covariant, is an application 

T: X • • • X ^(^)* X ^(^) X • • • X — > R 

•V ^ "V" 

A: times i times 

(cji,...,^fc,-j;i,...,^£) I — \ T{ui,...,Uk,vi,...,ve) 

(1.10) 

that is linear with respect to each of its arguments. The integer k + i is called the valence of the 
tensor. Let us recall the canonical duality = which means that every vector 

V can be considered as a linear form on the space ^(^)*, via v : ^(^)* — )■ M, i— )■ {u^v). 
Accordingly a vector is a tensor of type (1,0). A linear form is a tensor of type (0, 1). A tensor 
of type (0, 2) is called a bilinear form. It maps couples of vectors to real numbers, in a linear 
way for each vector. 

Given the natural basis {da) associated with some coordinates (x") and the corresponding 
dual basis (dx"), we can expand any tensor T of type (A;,£) as 

T = T"^-"'^^^^^^^^ 4i ® . . . ® 4, ^ dx^i dx^^ I, (1.11) 



where the tensor product dai (S> ■ ■ ■ da^ ^ dx^^ (X" . . . (g) dx^^ is the tensor of type (/c, i) for 
which the image of {uji, . . . , Uk,Vi, ■ ■ ■ , ve) as in (1.10) is the real number 

k i 

i=\ j=i 

Notice that all the products in the above formula are simply products in M. The 4*^+^ scalar 
coefficients (l-H) are called the components of the tensor T with respect 

the coordinates (x"). These components are unique and fully characterize the tensor T. 

Remark : The notation and uja already introduced for the components of a vector v or a hnear form 
Lj are of course the particular cases (fc, £) ~ (1, 0) and (fc, £) = (0, 1) of the general definition given 
above. 



1.1.4 Metric tensor 

As for special relativity, the fundamental structure given on the spacetime manifold ^ of general 
relativity is the metric tensor g. It is now a field on at each point p E g{p) is a bilinear 
form acting on vectors in the tangent space S/'pi^): 

g{p): ^p(^)x^p(^) M 

{u,v) I — > g{u,v). 

It is demanded that 

• gi is symmetric: g{u,v) = g{v,u)\ 



1.1 Geometrical framework 



11 



• gr is non- degenerate: any vector u that is orthogonal to all vectors (Vi> € ^p{^), g{u, v) 
0) is necessarily the null vector; 

• g has a signature (—,+,+,+): in any basis of where the components (ga^) form 
a diagonal matrix, there is necessarily one negative component {goo say) and three positive 
ones (gii, g22 and 533 say). 

Remark : An alternative convention for the signature of g exists in the literature (mostly in particle 
physics): (+,—,—,—). It is of course equivalent to the convention adopted here (via a mere 
replacement of g by g). The reader is warned that the two signatures may lead to different signs 
in some formulas (not in the values of measurable quantities !). 

The properties of being symmetric and non-degenerate are typical of a scalar product, 
thereby justifying the notation 



u v := g{u, v) = g^/s u"v'^ 



(1.13) 



that we shall employ throughout. The last equality in the above equation involves the components 
5a/3 of ff) which according to definition (1.11), are given by 

g = gap dx^" fS) dx^ . (1.14) 

The isotropic directions of g give the local light cones: a vector v € ^(^) satisfying 

V ■ V = is tangent to a light cone and called a null or lightlike vector. U v ■ v < 0, the vector 

V is said timelike and if -y • i; > 0, i; is said spacelike. 



1.1.5 Covariant derivative 

Given a vector field v on it is natural to consider the variation of iJ between two neighbouring 
points p and p' . But one faces the following problem: from the manifold structure alone, the 
variation of v cannot be defined by a formula analogous to that used for a scalar field [Eq. (1.8)], 
namely diJ = v{p') — v{p), because v{p') and v{p) belong to different vector spaces: 3pi{^) 
and £^p{^) respectively (cf. Fig. 1.1). Accordingly the subtraction v{p') — v{p) is not well 
defined ! The resolution of this problem is provided by introducing an affine connection on 
the manifold, i.e. of a operator V that, at any point p G associates to any infinitesimal 
displacement vector pf>' E ^(^) a vector of =5p(^) denoted ^p^''^ that one defines as the 
variation of v between p and p': 

dv = V^.v. (1.15) 

The complete definition of V allows for any vector u E =3^(^) instead of and includes all 
the properties of a derivative operator (Leibniz rule, etc.), which we shall not list here (see e.g. 
Wald's textbook [JO 1]). One says that v is transported parallelly to itself between p and p' 
iff dv = 0. 

The name connection stems from the fact that, by providing the definition of the variation dv, 
the operator V connects the adjacent tangent spaces ,'^p{^^) and ^p/{^). From the manifold 
structure alone, there exists an infinite number of possible connections and none is preferred. 
Taking account the metric tensor g changes the situation: there exists a unique connection, 
called the Levi-Civita connection, such that the tangent vectors to the geodesies with respect 
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to g are transported parallelly to themselves along the geodesies. In what follows, we will make 
use only of the Levi-Civita connection. 

Given a vector field v and a point p £ we can consider the type (1, 1) tensor at p denoted 
by 'Vv{p) and defined by 



By varying p we get a type (1,1) tensor field denoted Vi? and called the covariant derivative 
of V. 

The covariant derivative is extended to any tensor field by (i) demanding that for a scalar 
field V/ := d/ and (ii) using the Leibniz rule. As a result, the covariant derivative of a tensor 
field T of type {k, t) is a tensor field VT of type (fc, £ + 1). Its components with respect a given 
coordinate system (x") are denoted 

V7r"^-%,...ft:=(VT)— (1.17) 

(note the position of the index 7 !) and are given by 

rpai...aj, _ O rpai...ak i \^ rp<^i--- cr ...ak 

1=1 



where the coefficients T"^^ are the Christoffel symbols of the metric gi with respect to the 
coordinates (x"). They are expressible in terms of the partial derivatives of the components of 
the metric tensor, via 



where g"'^ stands for the components of the inverse matrix^ of {gap)'- 

g'^^g.p = (1.20) 

A distinctive feature of the Levi-Civita connection is that 

fl.21) 



Vg = 



Given a vector field u and a tensor field T of type (k, i), we define the covariant derivative 
of T along u as 

V^T ■.= VT{^^_^,u). (1.22) 

fc+£ slots 

Note that V^T is a tensor of the same type {k,l) as T and that its components are 

iyuTr-''%,,„p^ = u^y,T-'-'''p,...p, (1.23) 



Since g is non-degenerate, its matrix {gafi) is always invertible, whatever the coordinate system. 
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1.2 Einstein equation 

General relativity is ruled by the Einstein equation: 



„ 1 „ 8ttG^ 

R--Rg = ^T 



(1.24) 



where R is the Ricci tensor associated with the Levi-Civita connection V (to be defined below), 
R := tTgR = g^'^R^j^ is the trace of R with respect to g and T is the energy-momentum tensor of 
matter and electromagnetic field'^ (to be defined below). The constants G and c are respectively 
Newton's gravitational constant and the speed of light. In the rest of this lecture, we choose 
units such that 

G=l and c=l. (1.25) 

The Ricci tensor i? is a part of the Riemann tensor 7?., a type (1,3) tensor which 
describes the curvature of V, i.e. the tendency of two initially parallel geodesies to deviate 
from each other. Equivalently the Riemann tensor measures the lack of commutativity of two 
successive covariant derivatives: 

V^V,^;" - V.V^t;" = 7^"^^, (1.26) 
The Ricci tensor R is defined as the following trace of the Riemann tensor: 

It is expressible in terms of the derivatives of the components of the metric tensor with respect 
to some coordinate system (x") as 

where the T'^^^ have to be replaced by their expression (1.19). Although this not obvious on the 
above expression, the Ricci tensor is symmetric : Rj^a = Rap- 

Having defined the left-hand side of Einstein equation (1.24), let us focus on the right-hand 
side, namely the energy -momentum tensor T. To define it, let us introduce a generic observer 
O in spacetime. It is described by its worldline which is a timelike curve in and the 
future-directed unit vector u tangent to (the so-called observer ^-velocity): 

uu = -l. (1.29) 

An important operator is the projector _L onto the 3-dimensional vector space orthogonal to 
u, which can be considered as the local rest space of the observer O. _L is a type (1, 1) tensor, 
the components of which are 

l-''p = d''^ + u''u(s. (1.30) 

In particular, _L"^ti'^ = u"' + u'^u/^u^ = u°' + n°( — 1) = 0. The energy-momentum tensor T is 
the bilinear form defined by the following properties which must be valid for any observer O: 



^Or any kind of field which are present in spacetime. 
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• the energy density measured by O is 

E = T{u, u) = T^y uf'u'' ; (1.31) 

• the momentum density measured by O is the 1-form p the components of which are 

p„ = -T^,7x^±^^ ; (1.32) 

• the energy flux "vector"^ measured by O is the 1-form (fi the components of which are 

^a = -T^uU>'l-\; (1.33) 

• the stress tensor measured by O is the tensor S the components of which are 

S^^ = T^.L^^^L'^p. (1.34) 

A basic property of the energy-momentum tensor is to be symmetric. Anyway, it could not be 
otherwise since the left-hand side of Einstein equation (1.24) is symmetric. From Eqs. (1.32) and 
(1.33), the symmetry of T implies the equality of the momentum density and the energy flux 
form^: p = ip . This reflects the equivalence of mass and energy in relativity. 

Another property of T follows from the following geometric property, known as contracted 
Bianchi identity: 

1 



V^(^i?aM- 2^5aA.J =0, (1.35) 

where := g^'^V^- The Einstein equation (1.24) implies then the vanishing of the divergence 
of T: ^ ^ 

(1.36) 



V^T,^ = 



This equation is a local expression of the conservation of energy and momentum. 

An astrophysically very important case of energy-momentum tensor is that of a perfect 
fluid: 



[e +p)u®u + pg ■, (1-37) 



where u is the 1-form associated by metric duality^ to a unit timelike vector field u representing 
the fluid 4-velocity, e and p are two scalar fields, representing respectively the energy density and 
the pressure, both in the fluid frame. Indeed, plugging Eq. (1.37) into Eqs. (1.31)-(1.34) leads 
respectively to the following measured quantities by an observer comoving with the fluid: 

E = £, Pa = 0, ,fa=0, Sai3 = P-i-al3- 

1.3 3+1 formalism 

The 3+1 decomposition of Einstein equation is the standard formulation used in numerical 
relativity [4, 15, 46]. It is also very useful in the study of rotating stars. Therefore we introduce 
it briefly here (see e.g. [46] for an extended presentation). 



4 



For an electromagnetic field, this is the Poynting vector. 
^Restoring units where c 7^ 1, this relation must read p — c~^ip. 
®The components of u are Ua ~ Qapu^ ■ 
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t+dt 



Figure 1.2: Foliation of the spacetime ^ by a family of spacelike hypersurfaces (I]t)igK. 
1.3.1 Foliation of spacetime 

The 3+1 formalism relies on a slicing of spacetime by a family of spacelike hypersurfaces {I]t)t&R- 
Hypersurface means that is a 3-dimensional submanifold of spacelike means that every 
vector tangent to Tit is spacelike and slicing means (cf. Fig. 1.2) 



with Tit C^Tf =0 for t ^ t' . Not all spacetimes allow for a global foliation as in Eq. (1.38), but 
only those belonging to the class of the so-called globally hyperbolic spacetimes^ . However 
this class is large enough to encompass spacetimes generated by rotating stars. 

At this stage, i is a real parameter labelling the hypersurfaces Tt, to be identified later with 
some "coordinate time". We may consider t as a scalar field on It gradient Vt = d/ is then 
a 1-form that satisfies (Vt, v) = for any vector v tangent to Tt- This is direct consequence of t 
being constant on Tt- Equivalently the vector Vt associated to the 1-form Vt by metric duality 
(i.e. the vector whose components are V^t = g°'^V^t) is normal to Tt- Vt ■ v = for any vector 
V tangent to Sj. Since Tt is spacelike, it possesses at each point a unique unit timelike normal 
vector n, which is future-oriented (cf. Fig. 1.2): 




(1.38) 



n • n = — 1. 



(1.39) 



The two normal vectors Vt and n are necessarily colinear: 



n = -NVt . 



(1.40) 



^Cf. [3fi] or § 3.2.1 of [Ki] for the precise definition. 
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The proportionality coefficient N is called the lapse function. The minus sign is chosen so that 
> if the scalar field t is increasing towards the future. 

Since each hypersurface is assumed to be spacelike, the metric 7 induced*^ by g onto 
is definite positive, i.e. 

yve^^pi^t), v^O^-f{v,v)>0. (1.41) 
Considered as a tensor field on the components of -f are given in terms of the components 
of the normal via 

7a/3 = 9al3 + ?^Q?^/3■ (1-42) 
Note that if we raise the ffi'st index via g"'^, we get the components of the orthogonal projector 
onto Sj [compare with (1.30)] : 

7";3 = 5"^ + n%^, (1.43) 

which we will denote by 7. 

Let D be the Levi-Civita connection in associated with the metric 7. It is expressible in 
terms of the spacetime connection V and the orthogonal projector 7 as 

I3l...l3q—I Ml ' P-Pl Pi ' Pq' P^<^^ Ul...Uq- {^■'^'^) 

1.3.2 Eulerian observer or ZAMO 

Since n is a unit timelike vector, we may consider the family of observers whose 4-velocity is ft. 
Their worldlines are then the field lines of n and are everywhere orthogonal to the hypersurfaces 
Sj. These observers are called the Eulerian observers. In the context of rotating stars or 
black holes, they are also called the locally non-rotating observers [ ] or zero-angular- 
momentum observers (ZAMO) [102] . Thanks to relation (1.40), the proper time r of a 
ZAMO is related to t by 

dT = Ndt, (1.45) 

hence the name lapse function given to N. The 4-acceleration of a ZAMO is the vector a := V^n. 
This vector is tangent to and is related to the spatial gradient of the lapse function by (see 
e.g. § 3.3.3 of [ ]( )] for details) 

a = ^lnAf = 7(VlnA^). (1.46) 

In general a 7^ 0, so that the ZAMO's are accelerated observers (they "feel" the gravitational 
"force"). On the other side, they are non-rotating: being orthogonal to the hypersurfaces Sj, 
their worldlines form a congruence without any twist. Physically this means that the ZAMO's 
do not "feel" any centrifugal nor Coriolis force. 

Let us express the energy-momentum tensor T in terms of the energy E density, the momen- 
tum density p and the stress tensor S, these three quantities being measured by the ZAMO: 
setting u = a and _L = 7 in Eqs. (1.31), (1-32) and (1.34), we get 

E = T^^n^n" (1.47) 
Pa = -r^.n'^7^„ (1-48) 

Sap = T^.raYp- (1-49) 

The associated 3+1 decomposition of the energy-momentum tensor is then: 

Tap = Enani3+pani3+naPp + Sap. (1.50) 



7 is nothing but the restriction of g to Et 
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t+dt 



Figure 1.3: Coordinates (x') on the hypersurfaces Sj: each line = const cuts across the foliation 
(Et)tgR and defines the time vector dt and the shift vector /3 of the spacetime coordinate system {x") = 



1.3.3 Adapted coordinates and shift vector 

A coordinate system (x") is said to be adapted to the foliation [T,t)teR iff = t. Then the 
triplet^ (x*) = {x^,x'^,x^) constitutes a coordinate system on each hypersurface S^, which we 
may cah spatial coordinates. Given such a coordinate system, we may decompose the natural 
basis vector dt into a part along n and a part tangent to S^: 



dt = Nn + l3 



with n-/3 = 0. (1.51) 



The spacelike vector /3 is called the shift vector. Indeed it measures the shift of the lines of 
constant spatial coordinates with respect to the normal to the hypersurfaces T,t (cf. Fig. 1.3). The 
fact that the coefficient of n in Eq. (1.51) is the lapse function N is an immediate consequence 
of {dt,dt) = 1 [Eq. (1.4)] and relation (1.40). As for any vector tangent to T^, the shift vector 
has no component along dt : 

0=(3'di. (1.52) 
Equation (1.51) leads then to the following expression for the components of n: 

„ /I /52 /?3\ 

The covariant components (i.e. the components of the 1-form associated with n by metric 
duality) are given by Eq. (1.40) and (1.6): 

n„ = (-iV, 0,0,0). (1.54) 

Remark : It is immediate to check on (1.53)-(1.54) that n is a timelike unit vector: n n = n^n^^ = 
{-N){l/N) = -l. 



Latin indices ...) run in {1,2,3}, whereas Greek indices (a,/3, ...) run in {0, 1,2,3}. 



18 General relativity in brief 



The components {ga/s) of the spacetime metric are expressible in terms of the components 
i'Jij) of the induced metric in Et, the components of the shift vector and the lapse function: 



gaB dx^ dx 



/3 



-N^dt^ + 7ij(dx* + I3'dt){dx^ + P^dt) 



:i.55) 



1.3.4 Extrinsic curvature 

The intrinsic curvature of the hypersurface T,t equipped with the induced metric 7 is given 
by the Riemann tensor of the Levi-Civita connection D associated with 7 (cf. § 1.2). On the 
other side, the extrinsic curvature describes the way is embedded into the spacetime {^,g). 
It is measurable by the variation of the normal unit vector n as one moves on Tit- More precisely, 
the extrinsic curvature tensor K is the bilinear form defined on by 

V(m, v) G X K{u, v) := -u ■ n. (1.56) 

It can be shown (see e.g. § 2.3.4 of [Ki]) that, as a consequence of ft being hypersurface- 
orthogonal, the bilinear form K is symmetric (this result is known as the Weingarten prop- 
erty). 

The components of K with respect to the coordinates (x*) in E^ are expressible in terms of 
the time derivative of the induced metric 7 according to 

( t - ^f-") -^ir"^-^^- - ■ <-^) 

stands for the components of the Lie derivative of 7 along the vector field j3 (cf. Ap- 
pendix A) and the second equality results from the 3-dimensional version of Eq. (A. 8). The 
trace of K with respect to the metric 7 is connected to the covariant divergence of the unit 
normal to Et : 

K := Y^Kij = -V^. (1.58) 



1.3.5 3+1 Einstein equations 

Projecting the Einstein equation (1.24) (i) twice onto E^, (ii) twice along n and (iii) once on Ef 
and once along n, one gets respectively the following equations [4, 15, 46] : 



dt 



C^Kij = -D.DjN + N {^Rij + KKij - 2KikK^j + 47r [{S - E)j,j - 25^,]} 



:i.59) 
:i.6o) 

:i.6i) 



In this system, pi and Sij are the matter quantities relative to the ZAMO and are defined 
respectively by Eqs. (1.47)-(1.49). The scalar S is the trace of S with respect to the metric 7: 
S = Sij. The covariant derivatives Di can be expressed in terms of partial derivatives with 
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respect to the spatial coordinates (x*) by means of the Christoffel symbols ^r*^-^, of D associated 
with (x*): 

n,D,N=^-^T\,f,, (1.62) 

D,K\ = ^ + ''r^,,K\ - ''V\,K\, (1.63) 
dK 

D.K = -, (1.64) 

The C0Kij are the components of the Lie derivative of the tensor K along the vector /3 (cf. 
Appendix A); according to the 3-dimensional version of formula (A. 8), they can be expressed in 
terms of partial derivatives with respect to the spatial coordinates (x*) : 

The Ricci tensor and scalar curvature of 7 are expressible according to the 3-dimensional analog 
of formula (1.28): 



»j _ " ^ ik _ 

Q^k Q^j ^ ^ V ^ kl ^ ik ^ Ij 

'R = Y''^Rv- (1-67) 



- i^u- + r »j r - r r ; . (1.66) 



Finally, let us recall that Kij is related to the derivatives of ■~fij by Eq. (1.57). 

Equations (1.60) and (1.61), which do not contain any second order time derivative of jij, 
are called respectively the Hamiltonian constraint and the momentum constraint. 

The 3+1 Einstein equations are at the basis of the formulation of general relativity as a 
Cauchy problem and are much employed in numerical relativity (for reviews see e.g. [36, 46, 4]). 
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Chapter 2 

Stationary and axisymmetric 
spacetimes 
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Having reviewed general relativity in Chap. 1, we focus now on spacetimes that possess two 
symmetries: stationarity and axisymmetry, in view of considering rotating stars in Chap. 3. 

2.1 Stationary and axisymmetric spacetimes 

2.1.1 Definitions 
Group action on spacetime 

Symmetries of spacetime are described in a coordinate-independent way by means of a (sym- 
metry) group acting on the spacetime manifold ^ . Through this action, each transformation 
belonging to the group displaces points within ^ and one demands that the metric g is invariant 
under such displacement. More precisely, given a group G, a group action of G on ^ is an 
application^ 

$ : G X ^ — > ^ 

{9:P) ' — > ^{g,p) =■■ 9{P) 

such that 

• V((7, h) G G^, \/p G g{h{p)) = gh{p), where gh denotes the product of g hy h according 
to the group law of G (cf. Fig. 2.1); 

• if e is the identity element of G, then Vp G e{p) = p . 

The orbit of a point p G ^ is the set {g{p), g G G} C i.e. the set of points which are 
connected to p by some group transformation. 

^Do no confuse the generic element g of the group G with the metric tensor g. 
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P' 





Figure 2.1: Group action on the spacetime manifold 



An important class of group actions are those for which G is a one-dimensional Lie group 
(i.e. a "continuous" group). Then around e, the elements of G can be labelled by a parameter 
t G M, such that gt=o = e. The orbit of a given point p £ ^ under the group action is then 
either {p} (case where p is fixed point of the group action) or a one-dimensional curve of ^ . In 
the latter case, t is then a natural parameter along the curve (cf. Fig. 2.2). The tangent vector 
corresponding to that parameter is called the generator of the symmetry group associated 
with the t parametrization. It is given by 



where dx is the infinitesimal vector connecting the point p to the point gdt{p) (cf- § 1.1.1 and 
Fig. 2.2). The action of G on ^ in any infinitesimal neighbourhood of p amounts then to 
translations along the infinitesimal vector dt$,. 

Stationarity 

A spacetime {^,g) is said to be stationary iff there exists a group action on ^ with the 
following properties: 

1. the group G is isomorphic to (R, -|-), i.e. the group of unidimensional translations; 

2. the orbits are timelike curves in 

3. the metric is invariant under the group action, which is translated by 



where ^ is the generator of G associated with some parameter t of G and g denotes 

the Lie derivative of g along the vector field The Lie derivative measures the variation 
of g along the field lines of ^, i.e. the variation under the group action, and is defined in 
Appendix A. 




(2.3) 
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Figure 2.2: Orbit of a point p under tlie action of a one-dimensional Lie group, the elements gt of which 
being labelled by the parameter ^ e M. The vector ^ = dx/dt is the group generator associated with this 
parameter. 



Remark 1: The property (2.3) expresses the fact that if two vectors u and v are invariant under trans- 
port along the field lines of ^, then their scalar product u ■ v is also invariant. Indeed, thanks to 
the Leibniz rule: 

V^{u-v) = C^{u-v) = C^[g{u,v)] = {C^g){u,v) + g{C^u,v) + g{u, C^v), 
so that if C^u ~ and C^v = hold, C^g = implies V^(tZ ■ v) ^ 0. 

Remark 2: The property (2.3) is independent of the choice of the generator ^. i.e. of the parametrization 
t of G. Indeed, under a change of parametrization 1 1— >■ t', inducing a change of generator ^ i— > 
the following scaling law holds: 

dt ^ 
C^,g = ^C^g. 

Expressing the Lie derivative g according to Eq. (A. 11) and using V^ga/s = [property (1.21)] 
as well as = Safi^^ shows that the condition (2.3) is equivalent to 



VaC/J + y^^a = . (2. 



This equation is known as Killing equation^. Accordingly, the symmetry generator ^ is called 
a Killing vector. 

For an asymptotically flat spacetime, as we are considering here, we can determine the Killing 
field uniquely ^ by demanding that far from the central object ^ • ^ — > —1. Consequently, the 
associated parameter t coincides with the proper time of the asymptotically inertial observer at 
rest with respect to the central source. In the following we shall employ only that Killing vector. 



^Named after the German mathematician Wilhelm Killing (1847-1923). 
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Staticity 



A spacetime g) is said to be static iff 



1. it is stationary; 



2. tlie Killing vector field ^ is orthogonal to a family of liypersurfaces. 



Remark : Broadly speaking, condition 1 means that nothing depends on time and condition 2 that 
there is "no motion" in spacetime. This will be made clear below for the specific case of rotating 
stars: we will show that condition 2 implies that the star is not rotating. 

Axisymmetry 

A spacetime g) is said to be axisymmetric iff there exists a group action on ^ with the 
following properties: 

1. the group G is isomorphic to SO (2), i.e. the group of rotations in the plane; 

2. the metric is invariant under the group action: 



X being the generator of G associated with some parameter (p of G and g denoting the 
Lie derivative of g along the vector field %; 

3. the set of fixed points under the action of G is a 2-dimensional surface of which we 
will denote by A. 

Carter [-32] has shown that for asymptotically fiat spacetimes, which are our main concern here, 
the first and second properties in the above definition imply the third one. Moreover, he has 
also shown that A is necessarily a timelike 2-surface, i.e. the metric induced on it by g has the 
signature (—,+). A is called the rotation axis. 

Remark : The 2-dimensional and timelike characters of the rotation axis can be understood by con- 
sidering A as the time development (history) of the "standard" one-dimensional rotation axis in a 
3-dimensional space. 

From the very definition of a generator of a symmetry group, the vector field x must vanish on 
the rotation axis (otherwise, the latter would not be a set of fixed points): 



C^g = 



(2.5) 



xIa = o- 



(2.6) 



In addition, as for ^, the condition (2.5) is equivalent to demanding that % be a Killing vector: 



(2.7) 



Given a S0(2) group action, we can determine uniquely the Killing vector x t>y demanding 
that the associated parameter ip takes its values in [0,27r[. In the following, we shall always 
employ that Killing vector. 
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Figure 2.3: Commutativity of the stationary and axisymmetric actions. 
2.1.2 Stationarity and axisymmetry 

We consider a spacetime g) that is both stationary and axisymmetric. Carter has shown in 
1970 [32] that no generality is lost by considering that the stationary and axisymmetric actions 
commute. In other words, the spacetime g) is invariant under the action of the Abelian group 
M X S0(2), and not only under the actions of M and S0(2) separately. Saying that the stationary 
and axisymmetric actions commute means that starting from any point p G moving to the 
point q = Tt{p) under the action of an element Tt of the group (M, +) and then displacing q 
via an element of SO (2) yields the same point as when performing the displacements in the 
reverse order (cf. Fig. 2.3): 

Vp G ^, Vt G M, V99 G [0, 27r[, R^{Tt{p)) = Tt{R^{p)). (2.8) 

This property can be translated in terms of the Killing vectors ^ and x associated respectively 
with the parameter t of (M, +) and the parameter (p of S0(2). Indeed (2.8) is equivalent to the 
vanishing of their commutator: 



[tx]=0, (2.9) 



the commutator being the vector whose components are 

[txr = - X^'lg = e'^V.x" - x^V.r. (2.10) 

An important consequence of the above commutation property is that the parameters t and 
ip labelling the two symmetry groups (M, +) and S0(2) can be chosen as coordinates on the 
spacetime manifold Indeed, having set the coordinates (t, ip) = (0, 0) to a given point p, 
(2.8) implies that we can unambiguously attribute the coordinates {t,ip) to the point obtained 
from p by a time translation of parameter t and a rotation of parameter cp, whatever the order 
of these two transformations (cf. Fig. 2.3). 

Let us complete {t,Lp) by two other coordinates {x^,x'^) = {r,9) to get a full coordinate 
system on ^: 



(a;") = (t,r,^,V,) . (2.11) 
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{r,9) are chosen so that the spatial coordinates (r,6,ip) are of spherical type: r G [0, +oo[, 
6 G [0, tt] and = or = vr on the rotation axis A (cf. Fig. 2.4). 

Remark : An alternative choice would have been {x^,x^) = {p,z) such that the spatial coordinates 
{p,z,ip) are of cylindrical type: p G [0, +oo[, z G M and p = on A. The relation between the two 
types of coordinates is p = rsinO and z ~ rcos9. 

By construction, the coordinate system (2.11) is such that the first and fourth vectors of the 
associated natural basis (cf. § 1.1.1) are the two Killing vectors: 



dt 



and 



(2.12) 



Consequently, in terms of tensor components in the coordinates (x"), the Lie derivatives with 
respect to ^ (resp. x) reduce to partial derivatives with respect to t (resp. ip) [cf. Eq. (A. 3)]. In 
particular, the stationary and axisymmetry conditions (2.3) and (2.5) are respectively equivalent 

to 



dg, 



'a(5 



dt 







and 



dg, 



al3 



dip 







(2.13) 



For this reason, t and ip are called ignorable coordinates. The four coordinates (x") are said 
to be adapted coordinates to the spacetime symmetries. 

There is some freedom in choosing the coordinates (x") with the above properties. Indeed, 
any change of coordinates of the form 



t' 



t + T{r,9) 
R{r, 6) 
e{r,9) 
ip + <l>ir,e) 



(2.14) 



where T, R, @ and ^> are arbitrary (smooth) functions'^ of {r,9), lead to another adapted coor- 



Of course, to preserve the spherical type of the spatial coordinates, the functions R and O have to fulfill 
certain properties. 
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Figure 2.5: Change of coordinate t' according to (2.14). The solid lines represent the hypersurfaces 
of constant t (St) and the dashed lines those of constant t' (SJ,). The dotted lines are the orbits of the 
stationarity action (generator ^). 

dinate system: 

dti = i and d^> = X- (2-15) 

It is important to realize that in the above equation, ^ and x are the same vectors than in 
Eq. (2.12). The change of coordinate t ^ t' according to (2.14) is merely a reparametrization of 
each orbit of the stationarity action — reparametrization which leaves invariant the associated 
tangent vector ^ (cf. Fig. 2.5). A similar thing can be said for ip i— )• ip' . 

Remark : As illustrated in Fig. 2.5, stationarity by itself does not introduce a privileged 3+1 slicing of 
spacetime (cf. § 1.3), for the change of t coordinate according to (2.14) changes the hypersurfaces 
of constant t. 

2.2 Circular stationary and axisymmetric spacetimes 
2.2.1 Orthogonal transitivity 

Regarding the components of the metric tensor, the properties (2.13) are a priori the only ones 
implied by the spacetime symmetries (stationarity and axisymmetry) . There is however a wide 
subclass of stationary and axisymmetric spacetimes in which, in addition to (2.13), we are allowed 
to set to zero five metric components gap: the so-called circular spacetimes. Furthermore, these 
spacetimes are much relevant for astrophysics. 

Let us first remark that having g^p = for some value of (a, (3) is an orthogonality condition: 
that of the vectors da and dp^ i.e. of the lines {x"' = const, a' ^ a} and {x^' = const, /?' 7^ 
Next, we note that, in the present context, there are privileged 2-dimensional surfaces in 
spacetime: the orbits of the ]RxS0(2) group action. They are called the surfaces of transitivity 
and denoted S^Ai where A is some label. The two Killing vectors ^ and x are everywhere tangent 
to these surfaces and, except on the rotation axis, they form a basis of the tangent space to 
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at each point. Once adapted coordinates = (t, r, 9, y?) are chosen, the surfaces of transitivity 
can be labelled by the value of (r, 9) since both r and 9 are fixed on each of these surfaces: 

,9'A = yre = {p^^ / x\p) = r, x'^ip) = 6}. (2.16) 

Remark : r6 is a mere label for the surfaces S^re', the latter depend only on the symmetry group 
M X S0(2) and not on the choice of the {r,9) coordinates. A coordinate change {r,9) i— ?> {r',9') 
according to (2.14) will simply result in a relabelling of the surfaces =5^re- 

Given the family of 2-surfaces (J^^e)) one may ask if there exists another family of 2-surfaces, 
{Mb) say, which are everywhere orthogonal to ^r9- If this is the case, the M x S0(2) group 
action is said to be orthogonally transitive [31] and the spacetime {^,g) to be circular. 
One may choose the coordinates {r,9) to lie in Mb, i-e. the label B to be tip since both t and 
(p are then fixed on each of these surfaces: 

Mb = Mt^ = {v^J^ I x^{p) = t, x^{p) = (2.17) 

The 2-surfaces M b are called the meridional surfaces. In case of orthogonal transitivity, the 
following metric components are identically zero, refiecting the orthogonality between J^j-e and 

gtr = 0, gt9 = 0, g^r = 0, g^e = 0. (2.18) 
The following theorem states under which conditions this is guaranteed: 

Generalized Papapetrou theorem: a stationary (Killing vector ^) and axisym- 
metric (Killing vector x) spacetime ruled by the Einstein equation is circular iff the 
energy-momentum tensor T obeys to 

X^T^H^X^^ = 0, 

where the square brackets denote a full antisymmetrization. 

This theorem has been demonstrated in the case of vacuum solutions (T = 0) by Papapetrou 
(1966) [85] and extended to the non- vacuum case by Kundt &; Triimper (1966) [68] and Carter 
(1969) [31, 33] (see also § 7.1 of Wald's textbook [104], § 7.2.1 of Straumann's one [98], or § 19.2 
of [94]). Note that Eqs. (2.19)-(2.20) are equivalent to 

VeYect{^,x) (2.21) 
WeYect{i,x), (2.22) 

where V (resp. W) is the vector whose components are V" = ^^T^^ (resp. W"' = x^T^") and 

Vect(^,x) denotes the vector plane generated by ^ and X- 

In the important case of a perfect fluid source, the energy-momentum tensor takes the form 
(1.37), so that 

er^H^X^^ = ie+p)en,u^''^^x^^ +PeS,H^X^^ 

= (e+p)e'^n^^[-eV+p|[^|V^ 



= (e+p)e%u["e^x^^. 



(2.19) 
(2.20) 
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Similarly, 



Since ^ and u are both timelike, we have ^^'n^ ^ 0. The circularity conditions (2.19)-(2.20) are 
therefore equivalent to 

y[a^/3^7] = 0, (2.23) 



i.e. to 



Vect(C,x) 



(2.24) 



Taking into account that = dt and X = d^p, the above condition is equivalent to li'' = and 
= 0, or 



with 



u 



It* dt 



(2.25) 



(2.26) 



The 4-velocity (2.25) describes a pure circular motion of the fluid around the rotation axis, 
hence the qualifier circular given to spacetimes that obeys the orthogonal transitivity property. 
In such case, there is no fluid motion in the meridional surfaces, i.e. no convection. Stationary 
and axisymmetric spacetimes that are not circular have been studied in Refs. [47, 18]. 



In all what follows, we limit ourselves to circular spacetimes. 



2.2.2 Quasi-isotropic coordinates 

In a stationary and axisymmetric spacetime {^,g) that is circular, we may use adapted coor- 
dinates (x") = {t,r,6,ip) such that {r,6) span the 2-surfaces Aitip which are orthogonal to the 
surfaces of transitivity S^re, leading to the vanishing of the metric components listed in (2.18). 
Moreover, we can always choose the coordinates (r, 9) in each 2-surface A4t<p such that the metric 
induced by g takes the form^ 

gab dx^dx'' = A^{dr'^ + r'^dO'^), (2.27) 

where A = A{r,9). Indeed, dr^ + r'^dO'^ is the line element of the 2-dimensional flat metric in 
polar coordinates and, in dimension 2, all the metrics are conformally related^ meaning that 
they differ only by a scalar factor A'^ as in (2.27). 

Remark : This is no longer true in dimension 3 or higher: in general one cannot write the metric line 
element as a conformal factor times the flat one by a mere choice of coordinates. 

Note that the choice of coordinates (r, 0) leading to (2.27) is equivalent to the two conditions: 

Qre = and gee = r'^grr- (2.28) 

The coordinates (t, r, 6, ip) with the above choice for (r, 0) are called quasi-isotropic coordi- 
nates (QI). The related cylindrical coordinates {t,p,z,(p), with p := rsin^ and z := rcos9, 

^Indices a and b take their values in {1, 2}. 
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are called Lewis- Papapetrou coordinates, from the work of Lewis (1932) [i 2] and Papapetrou 
(1966) [85]. 

Let us define the scalar function lo = uj{r,6) as minus the scalar product of the two Killing 
vectors ^ and normalized by the scalar square of X- 



x-x 



(2.29) 



As we will see below, the minus sign ensures that for a rotating star, w > (cf. Fig. 3.3). Since 
9tip=i-X and g^^ = X- X, we may write 

gtip = -^gifiip- (2.30) 

Besides, let us introduce the following function of {r,9): 

,2 9<fi<p 



y.2 gjj^2 I 



(2.31) 



Collecting relations (2.18), (2.27), (2.30) and (2.31), we may write the components of the metric 
tensor in the form 



dx"" dxl^ = -N'^df + A^{dr'^ + r'^dO'^) + sin^ e{dip - ujdtf 



where A^, A, B and uj are four functions of {r,9): 



N = iV(r, 9) 



1 


A = A{r,9) 


•> 


B = B{r,9) 




UJ = [x?(r, 6) 



Equivalently, in matrix form: 



\ 



N^ + B^u^^sin^B 








-ujB'^r'^ sin^ 



















^2^2 





—ujB'^r'^ sin^ 9 








sin2 9 



(2.32) 
(2.33) 

(2.34) 



The inverse of this matrix is 



0/3 





-Ar-2 








-io/N^ 


\ 
































v 


-00 /N^ 








(Srsin6')-2 -a;2/iV2 


/ 



(2.35) 



To prove it, it suffices to check that relation (1.20) holds. 

The local flatness of spacetime implies that the metric functions A and B must coincide on 
the rotation axis: 

A\^ = B\^. (2.36) 

Indeed, let us consider a small circle around the rotation axis at a fixed value of both t and 9, 
with 9 <^ 1, and centered on the point of coordinates (x") = (t,r, 0,0) (cf. Fig. 2.6). According 
to the line element (2.32), the metric length of its radius is i? ~ A{r,0)r9, whereas its metric 
circumference is C ~ B{r / cos 9,9) rsinB x 27r ~ 27rB{r,0) r9. The local flatness hypothesis 
implies that C = 2ttR. Would this relation not hold, a conical singularity would be present on 
the rotation axis. From the above values of C and R, we get A{r,0) = B{r,0), i.e. the property 
(2.36) for the half part of A corresponding to 9 = 0. The demonstration for the second part 
(9 = tt) is similar. 
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Figure 2.6: Radius R and circumference C of a small circle around the rotation axis {6 ^ 1). 



/3* = (0,0, -a;) 



and (iii) the following expression of the induced metric in the hypersurfaces : 



2.2.3 Link with the 3+1 formahsm 

In terms of the 3+1 formalism introduced in § 1.3, the comparison of (2.32) with (1.55) leads 
immediately to (i) the identification of as the lapse function, (ii) the following components of 
the shift vector: 

(2.37) 

(2.38) 
(2.39) 

(2.40) 



■yij dx'dx^ = A\dr^ + r^dO^) + B^r^ sin^ 9 dip 



Therefore 
and 



-iij = diag [A^, A^r^, B^r^ shi^ 6) 



7- 



diag ( 



1 



A2' ^2^2' B^r^sin^e 



Remark : Comparing with (2.34)-(2.35), note that = 7^- but ^ 7*-' . 

Moreover, since in the present case dt = ^ and (3 = —ujO^p = —OJX, relation (1.51) becomes 



^ = Nfi — LUX 



(2.41) 



In particular, if t^; = 0, then ^ is colinear to n and hence orthogonal to the hypersurfaces E^. 
According to the definition given in § 2.1.1, this means that for u = 0, the spacetime is static. 
The reverse is true, assuming that in case of staticity, the axisymmetric action takes place in the 
surfaces orthogonal to ^. We may then state 



uj = 



g) static 



(2.42) 
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Let us evaluate the extrinsic curvature tensor K of the hypersurfaces T,t (cf. § 1.3.4). Since 
d-fij/dt = and = except for = -w, Eq. (1.57) leads to 



1 
2N 



CO 







duj \ 



Since (7jj) is diagonal [Eq. (2.39)] and duj/d(p = 0, we conclude that all the components of K 
vanish, except for 



K 


— K — 






2N 


dr 






= KpQ = 


B'^r'^ svv? i 


3doj 


2N 


de 



In particular, the trace of K defined by (1.58) vanishes identically: 

1 







K 

2 a 



(2.43) 
(2.44) 

(2.45) 



This is equivalent to the vanishing of the divergence of the unit normal n to [cf. Eq. (1.58)]. 

is then of maximal volume with respect to nearby hypersurfaces (see e.g. § 9.2.2 of Ref. [46] for 
more details). This is fully analogous to the concept of minimal surfaces in a 3-dimensional Eu- 
clidean space (the change from minimal to maximal being due to the change of metric signature) . 
For this reason, the foliation (Ej) with K = is called a maximal slicing. 

The quadratic term KijK^^ which appears in the Hamiltonian constraint (1.60) is 



= 2f^^^^{Kr^f + 2^''^^^{Ke^)\ 
where to get the last line we used the fact that 7*-' is diagonal. We thus get 

sin^ e 



duduj, 



2A^N^ 

with the following short-hand notation for any scalar fields u and v : 



^ g dudv 1 du dv 
dr dr dO 89 



(2.46) 



(2.47) 
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3.1 General framework 

We focus at present on the case of a single rotating star in equilibrium. The corresponding 
spacetime {^,g) is then stationary and asymptotically flat. In addition, it is reasonable to 
assume that the spacetime is axisymmetric. This has been shown to be necessary for viscous 
and heat-conducting fluids by Lindblom (1976) [78]. For perfect fluids, the general argument is 
that a rotating non-axisymmetric body will emit gravitational radiation, and therefore cannot 
be stationary. In this respect the situation is different from the Newtonian one, where stationary 
non-axisymmetric rotating configurations do exist, the best known example being that of the 
Jacobi ellipsoids. 

In addition of being axisymmetric, we will assume that the spacetime is circular (cf. § 2.2.1), 
since this is relevant for a perfect fluid rotating about the axis of symmetry [cf. Eq. (2.25)]. 
We may then use the QI coordinates (x") = (t,r,9,ip) introduced in § 2.2.2. Consequently, the 
metric tensor is fully specified by the four functions A^, A, B and a; listed in (2.33). In this 
chapter we derive the partial differential equations that these functions have to fulfill in order 
for the metric to obey the Einstein equation (1.24), as well as the equations of equilibrium for a 
perfect fluid. We shall also discuss the numerical method of resolution. 
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3.2 Einstein equations in QI coordinates 
3.2.1 Derivation 

We shall derive the equations within the framework of the 3+1 formalism presented in § 1.3. 
Let us first consider the trace of the 3+1 Einstein equation (1.59) in which we make use of the 
Hamiltonian constraint (1.60) to replace + by IBttE + KijK^^ , to get 



8K 8K 

-Q^-P'7^ = -D.D^N + N [4niE + S) + K,,K^^] . 



(3.1) 



In the present context, K = and K^K^^ is given by Eq. (2.46). The Laplacian of N can be 
evaluated via the standard formula: 



(3.2) 



where ^77 is the square root of the determinant of the components ■jij of the metric 7. From 
expression (2.39), we get 

^ = A'^Br'^ sine. (3.3) 
Taking into account expression (2.40) for 7*-', we conclude that (3.1) is equivalent to 

1 ON 



92 AT 2 ON 1 AT 



2 56*2 r2tan6' 89 



A^A^N{E + S) + ^^^^^^^^8u}8oj -din B8N. (3.4) 



Let us now consider the Hamiltonian constraint (1.60). On the right-hand side, we have 
= and KijK^^ given by Eq. (2.46). There remains to compute the Ricci scalar ^R. This 
can be done by means of formulas (1.66)-(1.67). One gets 



'R 



8^ A 18A 1 (dAy ^ 1 8^ A 



Qy.2 r 8r A \ 8r 



1 



8A\- 



1 f8'^B Z8B 1 8'^B 



r2 9^2 ^2^ yQQ ) 

2 8B^ 



B V 9r2 r 8r 



2 86"^ r2 tan 6 89 



(3.5) 



The Hamiltonian constraint (1.60) thus becomes 

2 1 q2 



1 



92 A \^_}_(^y 1 8'^A 1 f8A\' 



B'^r^ sin2 9 ^ ^ 
H 8uj8uj 



4Ar2 



dB 



1 fd'^B 3 8B 1 8'^B 2 
B \ 8r^ r 8r 89"^ tan 9 89 



-SttA^E. 



(3.6) 



Moving to the momentum constraint (1.61), we first notice that DiK = 0. We may use the 
standard formula to evaluate the divergence of a symmetric tensor (such as K): 



1 8 

8x^ 



2 8x^ 



(3.7) 
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In the present case, the last term always vanishes since K^^ is non-zero only for non-diagonal 
terms (jk=r(p, ipr, 6(p or (p9) and 'jji^ is diagonal [Eq. (2.39)]. Since in addition d/d(p = 0, we 
are left with 



(3.8) 



Now, K^j. = K^j. = and K''^ = K^g = 0. We thus conclude that the first two components of 
the momentum constraint (1.61) reduce to 

Pr = and pq = 0. (3-9) 

For the third component, we use Eq. (3.8) for i = (p, taking into account that K^^ = j^^Kj^ = 
Y'^Kr^ = A-'^Kr^ and K^^ = Kj^ = -/^^Kq^ = A-'^r-^Ke^ with the values (2.43)-(2.44) for 
Kr:^ and Kg^. We thus obtain 



N 

S3 



sin0 d 



+ 



d 



r sm 



do 



sm 



-167r 



NA^ p^ 



r sm ( 



(3.10) 



Finally, let us consider the ipip component of the 3+1 Einstein equation (1.59). We have 
K^ip = and, according to Eq. (1.65), 







0. 



o^p dip 



Thus the left-hand side of the ipp component of Eq. (1.59) vanishes identically. On the right-hand 
side, we have DjN = dN/dx^ , 

ON 3 dN ON 



D D N = — — -3r^ 

dtp dip dx^ 



dr 







2 • 2 

r sm ( 



^2 

A^ 



1 dB l\ dN 1 



B dr 



dr 



2 • 2 , 

-r sm I 



B /d^B 3dB 1 d^B 



AP- V dr'^ ^ r dr 



and 



1 

l2 



(K ^2 



+ 



+ 



I dB 1 \ dN 
B~de ^ tan^y W 

2 dB' 



+ 



2 56*2 r2tan6' dO 



1 



B^r^ sin^ 9 
;{K^ef = — , .n,.n dujduj. 



LP ^^Lpr-- LP I -^Lpu^^ LP ^2V"'-</"''' ' y^2^2 ^""'/"'^ 4^2 iV2 

In addition, the matter term is 

{S - ^)7w - = iS\ + S% + S''^ - Eh^^ - 2-i^^S\ = r2 sin2 6 B\S\ + S\ - 5'^^ - E). 

Accordingly, after division by r"^ siv? 9 NB"^ /A^ , the (ptp component of the 3+1 Einstein equation 
(1.59) writes 



dB 



1 (d'^B 3dB 1 d'^B 2 
B \ dr"^ r dr r-^ d9'^ r-^ tan 9 d9 



+ ^dBdN + ^^^^d.du 



BN 



2iV2 



+ 



1 fdN 1 dN\ _ .,2(^r ,c.e 



E 



(3.11) 
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At this stage, we have four equations, (3.4), (3.6) (3.10) and (3.11), for the four unknowns 
N, A, B and u. We construct hnear combinations of these equations in which classical elliptic 
operators appear. First, multiplying Eq. (3.4) by Eq. (3.11) by NB and adding the two yields 



3 d 



1 a2 



{NB) + 



d 



tan 9 89 



(NB) = MNA^B{S\ + 5^ 



Second, dividing Eq. (3.4) by A^, adding Eq. (3.6) and subtracting Eq. (3.11) yields 
1 d'^N 1 dN 1 d'^N 1 d'^A 1 dA 1 d'^A 



1 dN 



1 d'^N 1 d'^A 

A dr"^ rA dr r'^A 89"^ 



SttA^S'p^ + —.dAdA + —5 du:duj. 



Reorganizing slightly Eqs. (3.4), (3.10), (3.12) and (3.13), we obtain the final system: 



A3U = A^A\E + S)+ 



dioduj — dud{v + In B) 



^^.■^{ujr sin ( 



-167r- 



NA^ 



Pip 



B^ 



r sm ( 



+ r sin 9 dujd{u — 3 In i?) 



A2 [{NB - l)rsin6l] = STiNA^Brsm9{S\. + S\) 



/S.2{\nA + v) 



^ttA^S'p^ + 



SB^r^ sin^ 
4iV2 



■ dudio — dvdu 



where the following abbreviations have been introduced: 



A2 
A3 
A3 



IniV 



1^ 1 92 

Qj.2 J. Qj. J.2 QQ2 

^2 2d_ 1 92 

Qj.2 J. Qj. J.2 QQ2 



+ 



1 



d 



' tan 9 89 



A3 



1 



r2 sin2 I 



(3.12) 



(3.13) 

(3.14) 

(3.15) 
(3.16) 
(3.17) 

(3.18) 
(3.19) 

(3.20) 
(3.21) 



Terms of the form 8vdv have been defined by (2.47). 

The operator A2 is nothing but the Laplacian in a 2-dimensional flat space spanned by the 
polar coordinates {r,9), whereas A3 is the Laplacian in a 3-dimensional flat space, taking into 
account the axisymmetry {8/d(p = 0). A deeper understanding of why these operators naturally 
occur in the problem is provided by the (2+l)+l formalism developed in Ref. [47]. 



3.2.2 Boundary conditions 

Equations (3.14)-(3.17) forms a system of elliptic partial differential equations. It must be sup- 
plemented by a set of boundary conditions. Those are provided by the asymptotic flatness as- 
sumption: for r —7- -|-oo, the metric tensor tends towards Minkowski metric 77, whose components 
in spherical coordinates are 



rjafi dx^ = -df + dr^ + r^d9^ + sin^ 9 d^"^. 



(3.22) 
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Comparing with (2.32), we get immediately the boundary conditions: 



N 
A 
B 



1 
1 
1 




when r 



-oo. 



(3.23) 



Remark : Spatial infinity is the only place where exact boundary conditions can be set, because the 
right-hand sides of Eqs. (3.14)-(3.17) have non-compact support, except for (3.16). This contrasts 
with the Newtonian case, where the basic equation is Poisson equation for tlie gravitational poten- 
tial $ [Eq. (3.43) below]: 

A3$ = 47rGp. (3.24) 

For a star, the right-hand side, involving the mass density p, has clearly compact support. The 
general solution outside the star is then known in advance: it is the (axisymmetric) harmonic 
function 

EPt (cos 9) 



$(r, ( 



(3.25) 



where Pi is the Legendre polynomial of degree i. In this case, one can set the boundary conditions 
for a finite value of r and perform some matching of $ and 5$ / dr to determine the coefficients ai . 

From the properties of the operators A3 and A3 and the boundary conditions (3.23), one can 
infer the following asymptotic behavior of the functions N and oj: 



N 



1 + ^ + 



1 



K 



1 



where Kq and Ki are two constants. In Chap. 4, we will show that 

Ko = -M and Ki = 2J, 



(3.26) 



(3.27) 



(3.28) 



M and J being respectively the mass of the star and its angular momentum [cf. Eqs. (4.16) and 
(4.41)]. 



3.2.3 Case of a perfect fluid 

Since we are interested in rotating stars, let us write the matter source terms in the system 
(3.14)-(3.17) for the case of a perfect fluid [cf. Eq. (1.37)]. Due to the circularity hypothesis, 
the fluid 4-velocity is necessarily of the form (2.25): u = m*(^ -|- 0,x)- We may use Eq. (2.41) to 
express ^ in terms of the unit timelike vector n (the 4-velocity of the ZAMO) and ojx. We thus 
get 

u=[Nu^) n+^{9.-uj)x . 



Now, according to Eq. (1.54), n • x 



n 



N 



(3.29) 



0. Therefore, Eq. (3.29) constitutes an 



orthogonal decomposition of the fluid 4-velocity with respect to the ZAMO 4-velocity n: 



u 



T{n + U 



(3.30) 
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with 



r = Nu' 



and 



u = -{n-uj)x 



(3.31) 



r is the Lorentz factor of the fluid with respect to the ZAMO and the vector U is the fluid 
velocity (3- velocity) with respect to the ZAMO (see e.g. Ref. [ ll] for details). U is a spacelike 
vector, tangent to the hypersurfaces Tit- Let us define 



U := — (17 — uj)r sin ( 



Since X' X = B^r^ sin^ 6, we have 



U -U = U\ 



(3.32) 
(3.33) 



From n • n = —1 and n • ?7 = 0, Eq. (3.30) yields the geometric interpretation of the Lorentz 
factor as (minus) the scalar product of the ZAMO 4-velocity with the fluid one: 



-n ■ u 



(3.34) 



The normalization condition u u = —1 expressed in terms of (3.30) leads to the standard relation 

= (l_[/2)-V2 I P 35^ 



The quantities E, p^p, S^^, S^g and S^,p which appear in the right-hand of Eqs. (3.14)-(3.17) 
are computed by means of Eqs. (1.47)-(1.49), using the form (1.37) of the energy momentum 
tensor. We get, using (3.30), 

E = T^uU^n" = {e+p) u^n^" u^n" +p g^^^n^n^ = T'^{e +p) -p, 



= -T^un^'-i^^ = -{e + p) u^jf^pUiyn^+p gf,i,n''-f''p = T'^{e + p)Up, 

ru^ 

5°^ = T^,7^-7'^/3 = {e+p)u^^>^»u,Yp+pg^.u7''''Y^ = r'{e + p)U''Uis+pjy 

Now, from Eq. (3.31), C/" = (0, 0, 0, C/"^), and from Eq. (3.33), U^Vf = U"^ and = Brsm9U. 
Besides 7*^^ = 7^^ = j^ip = 1. Accordingly, the above results can be expressed as 



E = T^{e + p)-p 



p^p = B{E + p)Ur sine 



S\ =p 


5 


S\=P 





S'^=p+{E + p)U' 



S = 3p+{E+p)U^ 



(3.36) 
(3.37) 
(3.38) 



From the above expressions, the term E + S in. the right-hand side of Eq. (3.14) can be written 

l + C/2 



E + S 



l-C/2 



{e+p) + 2p. 



(3.39) 
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3.2.4 Newtonian limit 

In the non-relativistic limit, v tends to the Newtonian gravitational potential <^ (divided by (?\. 



Newt. 



Moreover, Eq. (3.32) reduces to 



C7 = Or sin Q 



Newt. 



and [/^ ^ 1, p ^ e and e ~ pc^, where p is the mass density. Accordingly, (3.39) gives 



(3.40) 
(3.41) 
(3.42) 



Since in addition A ~ 1 and all the quadratic terms involving gradient of the metric potentials 
tend to zero, the first equation of the Einstein system (3.14)-(3.17) reduces to (after restoring 
the G and c's) 

' ' (3.43) 



A3<I> = 47rGp 



Newt. 



The other three equations become trivial, of the type "0 = 0". We conclude that, at the non- 
relativistic limit, the Einstein equations (3.14)-(3.17) reduce to the Poisson equation (3.43), which 
is the basic equation for Newtonian gravity. 



3.2.5 Historical note 

The Einstein equations for stationary axisymmetric bodies were written quite early. In 1917, 
Weyl [JU")] treated the special case of static bodies, i.e. with orthogonal Killing vectors ^ and x : 
or equivalently with cj = [cf. Eq. (2.42)]. The effect of rotation was first taken into account by 
Lense and Thirring in 1918 [71], who derived the first order correction induced by rotation in the 
metric outside a spherical body. The first full formulation of the problem of axisymmetric rotating 
matter was obtained by Lanczos in 1924 [(i!)] and van Stockum in 1937 [97]. They wrote the 
Einstein equations in Lewis-Papapetrou coordinates (cf. § 2.2.2) and obtained an exact solution 
describing dust (i.e. pressureless fluid) rotating rigidly about some axis. However this solution 
does not correspond to a star since the dust fills the entire spacetime, with an increasing density 
away from the rotation axis. In particular the solution is not asymptotically flat. Rotating stars 
have been first considered by Hartle and Sharp (1967) [GO], who formulated a variational principle 
for rotating fluid bodies in general relativity. Hartle (1967) [57] and Sedrakian & Chubaryan 
(1968) [90] wrote Einstein equations for a slowly rotating star, in coordinates different from the 
QI ones^. For rapidly rotating stars, a system of elliptic equations similar to (3.14)-(3.17) was 
presented (in integral form) by Bonazzola &; Machio (1971) [23]. It differs from (3.14)-(3.17) 
by the use of Ingu instead of IniV = i/. Bardeen & Wagoner (1971) [14] presented a system 
composed of three elliptic equations equivalent to Eqs. (3.14), (3.15) and (3.16), supplemented 
by a first order equation of the form 

^{lnA + u) = S{r,e), (3.44) 

where S{r,6) is a complicated expression containing first and second derivatives of i', A, B and 
cj and no matter term [see Eq. (4d) of Ref. [•■iO] (9 lines !) or Eq. (11.16) of [J I]]. Equation (3.44) 



^Namely coordinates for which goe — g^^/ sin^ 9 instead of (2.28). 
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results from the properties R^r = RbqIt^ and R^q = of the Ricci tensor. Note that these are 
verified only for an isotropic energy- momentum tensor, such as the perfect fluid one [Eq. (1.37)]. 
For anisotropic ones (like that of an electromagnetic field), Eq. (3.44) should be modified so as to 
include the anisotropic part of the energy-momentum tensor on its right-hand side. The full sys- 
tem of the four elliptic equations (3.14)-(3.17) has been first written by Bonazzola, Gourgoulhon, 
Salgado & Marck (1993)2 [22]. 



3.3 Spherical symmetry limit 
3.3.1 Taking the limit 

In spherical symmetry, the metric components g^fi simplify significantly. Indeed, we have seen 
that A and B must be equal on the rotation axis [Eq. (2.36)]. But, if the star is spherically 
symmetric, there is no privileged axis. Therefore, we can extend (2.36) to all space: 



A = B 



spher. sym. 



(3.45) 



Besides, to comply with spherical symmetry, the parameter in the decomposition (2.25) of the 
fluid 4-velocity must vanish. Otherwise, there would exist a privileged direction around which 
the fluid is rotating. Hence we have 



Q = 



spher. sym. 



(3.46) 



and u is colinear to the Killing vector ^. From Eqs. (3.37) and (3.32) with $7 = 0, we have 
Pip oc Lo. Equation (3.15) is then a linear equation in cj. It admits a unique solution which is 
well-behaved in all space: 

■ ■ (3.47) 



= 



spher. sym. 

According to (2.42), this implies that the spacetime is static. 

Properties (3.46) and (3.47), once inserted in Eq. (3.32) leads to the vanishing of the fluid 
velocity with respect to the ZAMO: 



[/ = 



spher. sym. 



Thanks to (3.45) and (3.47), the spacetime metric (2.32) becomes 



g^^ dx" dx^ = -N{rfdt^ + A{rf{dr'^ + r'^dO'^ + sin^ 9 d^'^ 



spher. sym. 



(3.48) 



(3.49) 



The coordinates (t, r, 6, if) are in this case truly isotropic, the spatial metric being conformal to 
a flat one. This justify the name quasi-isotropic given to these coordinates in the general case. 

Since U = implies F = 1, the Einstein equations (3.14)-(3.17) with a perfect fluid source 
and spherical symmetry reduce to 



d'^u ^2dv 4 _|_ 2 -J dv d 

dr'^ r dr dr dr 



{v + \nA) 



(3.50) 



spher. sym. 



In the oris inal article [22], the functions A = {ABf^^ and B = {B/A)^^'^ were employed instead of A and B. 
The equations with the present values of A and B are given in Ref. [!!)]. 
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ar^ r ar 



spher. sym. 



-^(In A + z.) + i-^(ln^ + z.) = S^A^p - f 
ar^ r ar \ar J 



(3.51) 



(3.52) 



spher. sym. 



We have omitted Eq. (3.15) since, as seen above, it is triviaUy satisfied in spherical symmetry. 
Besides, we have used the equivalent form (3.12) of Eq. (3.16). 

3.3.2 Link with the Tolman-Oppenheimer-Volkoff system 

The reader famihar with static models of neutron stars might have been surprised that the spheri- 
cal limit of the Einstein equations (3.14)-(3.17) does not give the familiar Tolman-Oppenheimer- 
Volkoff system (TOV). Indeed, the TOV system is (see e.g. [55, 58]): 



dm ^ _2 
— — = 47rr e 
ar 



du 
df 
dp 
df 



^= 1 



2m{r) 



m[r) 



+ 47rrpJ 



, ^ du 



(3.53) 
(3.54) 
(3.55) 



which is pretty different from the system (3.50)-(3.52). In particular, the latter is second order, 
whereas (3.53)-(3.55) is clearly first order. However, the reason of the discrepancy is already 
apparent in the notations used in (3.53)-(3.55) : the radial coordinate is not the same in the 
two system: in (3.50)-(3.52), it the isotropic coordinate r, whereas in the TOV system, it is the 
areal radius f. With this radial coordinate, the metric components have the form 



g^p dx" dx^ = -N{rfdf + ( 1 



2m(r) 



dr'^ + f'^d9'^ + r'^ siv? 6 d</ . 



(3.56) 



This form is different from the isotropic one, given by (3.49). Note however that the coordinates 
= (t, f, 0, if) differ from the isotropic coordinates = (t, r, 6, if) only in the choice of r. From 
(3.56), it is clear that the area of a 2-sphere defined by t = const and f = const is given by 



A = 47rf^ 



(3.57) 



hence the qualifier areal given to the coordinate r. On the contrary, we see from (3.49) that the 
area of a 2-sphere defined by t = const and r = const is 



2„2 



A = 47rA{ryr 



(3.58) 



so that there is no simple relation between the area and r. 

Outside the star, the solution is the Schwarzschild metric (see § 3.3.3) and the coordinates 
(x") are simply the standard Schwarzschild coordinates. In particular, outside the star m(f) = 
M = const, where M is the gravitational mass of the star and N{f)'^ = (1 — 2M/f), so that 
we recognize in (3.56) the standard expression of Schwarzschild solution. The expression of this 
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solution in terms of the isotropic coordinates used here is derived in § 3.3.3. The relation 
between the two radial coordinates is given by 

/ \ 2 

f = r ( 1 + — j (outside the star). (3.59) 

Far from the star, in the weak field region [Mjr <^ 1), we have of course f c^i r. 

There exist two coordinate systems that extend the TOV areal coordinates to the rotating 
case: 

• the coordinates employed by Hartle &: Thorne [57, 61] and Sedrakyan &: Chubaryan [90] to 
compute stellar models in the slow rotation approximation (metric expanded to the order 
^'); 

• the so-called radial gauge introduced by Bardeen & Piran [J 3]; however the entire hypersur- 
face cannot be covered in a regular way by these coordinates, as shown in Appendix A 
of Ref. [22]. 

3.3.3 Solution outside the star 

Outside the star, e = and p = 0. Equation (3.51) reduces then to 



— (NA) + - — (NA) - — — 
dr"^ r dr r'^ dr 



0, 



from which we deduce immediately that 



d 2a 
dr 



— (NA) = — , with a = const. 



The factor 2 is chosen for later convenience. The above equation is easily integrated, taken into 
account the boundary conditions (3.23): 

NA = 1--^. (3.60) 

Besides, subtracting Eq. (3.50) from Eq. (3.52), both with e = and p = 0, we get the following 
equation outside the star: 

d^ ^ , 1 d ^ , 1 du du d , , 

—rlnA+-—lnA-—- = ——lnA. 3.61 

ar^ r dr r dr dr dr 

Now Eq. (3.60) gives, since v = InN, 

dv d , , 2a 
■lnA + 



dr dr r{r'^ — a) 

Inserting this relation into Eq. (3.61), we get 



(PA 2(r2 - 2a) dA 2a , , , 

TT + ^Tl TTt 27-2 S^ = ^- 3-62 
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The reader can check easily that the unique solution of this linear ordinary differential equation 
which satisfies to the boundary condition (3.23) requires a > and is 




(3.63) 



Relation (3.60) yields then 




(3.64) 



The constant i/a is actually half the gravitational mass of the star. Indeed the asymptotic 
expansion for r — ?• oo of the above expression is 



r J \ r 
Comparing with Eqs. (3.26) and (3.28), we get 

V^=^. (3.65) 

Replacing a, A and N by the above values in (3.49), we obtain the following form of the metric 
outside the star: 

/l_M\^ / 

gc.^dx'' dx^ = -\—-^\ '^*^+(^l + ^J [dr"^ + r'^{de'^ + Wedip^)\. (3.66) 

We recognize the Schwarzschild metric expressed in isotropic coordinates. Actually, we have 
recovered the Birkhoff theorem^: With spherical symmetry, the only solution of Einstein 
equation outside the central body is the Schwarzschild solution. 



Remark : There is no equivalent of the Birkhoff theorem for axisymmetric rotating bodies: for black 
holes, the generalization of the Schwarzschild metric beyond spherical symmetry is the Kerr metric^ 
and the generic metric outside a rotating star is not the Kerr metric. Moreover, no fluid star has 
been found to be a source for the Kerr metric (see e.g. [<">()]). The only matter source for the Kerr 
metric found so far is made of two counterrotating thin disks of collisionless particles [17]. It has 
also been shown that uniformly rotating fluid bodies are finitely separated from Kerr solutions, 
except extreme Kerr [81]; this means that one cannot have a quasiequilibrium transition from a 
star to a black hole. 



3.4 Fluid motion 

Having discussed the gravitational field equations (Einstein equation), let us now turn to the 
equation governing the equilibrium of the fiuid, i.e. Eq. (1.36). 

^More precisely, we have established Birkhoff theorem under the hypothesis of stationarity; however this 
hypothesis is not necessary. 
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3.4.1 Equation of motion at zero temperature 

We consider a perfect fluid at zero temperature, which is a very good approximation for a 
neutron star, except immediately after its birth. The case of finite temperature has been treated 
by Goussard, Haensel & Zdunik (1997) [51]. The energy-momentum tensor has the form (1.37) 
and, thanks to the zero temperature hypothesis, the equation of state (EOS) can be written as 



e 
V 



e(nb) 



(3.67) 
(3.68) 



where nb is the baryon number density in the fluid frame. 

The equations of motion are the energy-momentum conservation law (1.36) : 



V^T,^ = 



and the baryon number conservation law: 



V;,(nbn'^) = 0. 



(3.69) 



(3.70) 



Inserting the perfect fluid form (1.37) into Eq. (3.69), expanding and projecting orthogonally to 
the fluid 4-velocity u [via the projector _L given by (1.30)], we get the relativistic Euler equation: 



(e + p)n^V^u„ + (5^„ + u^'ua)V^p = 0. 
Now the Gibbs-Duhem relation at zero temperature states that 

dp = nb d/i. 



(3.71) 



(3.72) 



where is the baryon chemical potential, /i := de/dn]^. Moreover, thanks to the first law of 
Thermodynamics at zero temperature (see e.g. Ref. [45] for details), /i is equal to the enthalpy 
per baryon h defined by 



h :-- 



e +p 
nb 



(3.73) 



Thus we may rewrite (3.72) as dp = ny, dh, hence 

VaP = n\^V ah. 
Accordingly, Eq. (3.71) becomes, after division by nb, 

hui'V^Ua + {d'^a + u''ua)V^h = 0, 
which can be written in the compact form 



u^'V^ihua) + Vah = 



(3.74) 



Thanks to the properties u^u^ = — 1 and u^'VaU^ = (the latter being a consequence of the 
former), Eq. (3.74) can be rewritten in the equivalent form 



[V^,{hUa)-VaihUf,)]=0 



(3.75) 



3.4 Fluid motion 



45 



This makes appear the antisymmetric bihnear form fi^^ := Va{hup) — Vaihujs), called the 
vorticity 2- form. The fluid equation of motion (3.75) has been popularized by Lichnerowicz 
[75, 76] and Carter [34], when treating relativistic hydrodynamics by means of Cartan's exterior 
calculus (see e.g. Ref. [45] for an introduction). 

Equation (3.74) [or (3.75)] is more useful than the original Euler equation (3.71). It leads 
easily to the relativistic generalization of the classical Bernoulli theorem (§ 3.4.2) and to a first 
integral of motion for a rotating star (§ 3.4.3). 



3.4.2 Bernoulli theorem 



Equation (3.74) is an identity between 1-forms. Applying these 1-forms to the stationarity Killing 
vector ^ (i.e. contracting (3.74) with ,^"), we get successively 

ru''v^(/iu,)+rv./i = o, 



w^V^(/in^D - hu^u^'Vf.C = 0, 



w^V^(/in^r) = 0. 

Note that ^"Vyh = stems from the stationarity of the fluid flow and u^vySJ = from 
Killing equation (2.4). We thus have 



V^{hu • ^) = 



(3.76) 



In other words, the scalar quantity hu ■ ^ is constant along any given fluid line. This results 
constitutes a relativistic generalization of the classical Bernoulli theorem. To show it, let us first 
recast Eq. (3.76) in an alternative form. Combining Eq. (3.30) with Eq. (2.41), we get (using 
n • n = —1, n • X = and ft ■ U = 0) 



hu-^ = hV{n + U) ■ {Nn - ux) = -hTN ( 1 + — x • U 



(3.77) 



Accordingly, by taking the logarithm of —hu ■ ^, we obtain that property (3.76) is equivalent to 

(3.78) 



OJ 



H + u + lnT + ln [1 + jy-^ ' ^) ~ const along a fluid line 



where we have introduced the log-enthalpy 



H ■.= ln{ 

mb 



(3.79) 



mb being the mean baryon mass : ?nb — 1.66 x 10 kg. Defining the fluid internal energy 
density by 

eint ■■= e- m^rib (3.80) 
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we have [cf. (3.73)] 



At the Newtonian hmit, 



H = \n{l + 



eint + P 
TTibnb 



(3.81) 



gjnt + P 



«1, 



and mb^-b — p (the mass density), so that H tends towards the (non-relativistic) specific enthalpy: 



H ~ 



eint + P 

p 



(3.82) 



Newt. 



Besides, thanks to (3.35) and U"^ <^ 1, we have InF ~ Moreover we have already seen 

that 1/ tends towards the Newtonian gravitational potential $ [Eq. (3.40)]. In addition, at the 
Newtonian limit, a; = 0. Consequently, the Newtonian limit of (3.78) is 



[/2 

H + '^ + — 
2 



const along a fluid line 



(3.83) 



Newt. 



which is nothing but the classical Bernoulli theorem. 

The general relativistic Bernoulli theorem (3.76) has been first derived by Lichnerowicz in 
1940 [73, 74]. The special relativistic version had been obtained previously by Synge in 1937 
[99]. 

Remark : That a symmetry (here stationarity, via the Killing vector ^) gives rise to a conserved quantity 
[Eq. (3.78)] is of course a manifestation of the Noether theorem. We will encounter another one in 
§ 4.6.1. 

To establish (3.76) we have used nothing but the fact that the flow obeys the symmetry and 
^ is a Killing vector. For a flow that is axisymmetric, we thus have the analogous property 



V^(/iM • X) = 



(3.84) 



The quantity hu ■ x is thus conserved along the fluid lines. Thanks to (3.30) and n ■ x = we 
may write it as 

hu-x = hTx-U = hTB'^r^sm^eU'P. (3.85) 



At the Newtonian limit, h 



1, r 



1, i? — 7- 1 and ~ so that the conserved quantity 



along each fluid line is the specific angular momentum about the rotation axis, sin^ 



Remark : It is well known that, along a timelike geodesic, the conserved quantities associated with the 
Killing vectors ^ and x ^ ' € ^nd it ■ x [^2, 104]. Here we have instead hu ■ ^ and hu ■ x, for 
the fluid lines are not geodesies, except when p = 0. In the latter case, h reduces to a constant 
(typically mb, cf. Eq. (3.73) and Eq. (3.80) with Eint = 0) and we recover the classical result. 
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3.4.3 First integral of motion 

In the derivation of the BernoulH theorem, we have assumed a general stationary fluid. Let us 
now focus to the case of a rotating star, namely of a circular fluid motion around the rotation 
axis^. The fluid 4-velocity takes then the form (2.25): 



u = u k 



(3.86) 



with 



(3.87) 



Note that if is a constant, fc is a Killing vector: it trivially satisfies Killing equation since ^ 
and X both do. u is then colinear to a Killing vector and one says that the star is undergoing 
rigid rotation. If is not constant, the star is said to be in differential rotation. 

To derive the first integral of motion, it is convenient to start from the fluid equation of 
motion in the Carter-Lichnerowicz form (3.75). Thanks to the symmetry of the Christoffel 
symbols (1.19) with respect to their last two indices, we may replace the covariant derivatives 
by partial ones in Eq. (3.75): 



d 



{huo 



d 



0. 



dxi^ dx°' 

Assuming a circular motion, i.e. H of the form (3.86), we have = (n*, 0, 0, ?i*ri), so that 



(3.88) 



d 



{hUa) = 



0. 



Hence Eq. (3.88) reduces to (after division by n* 



or equivalently 



_d_ 



[h{ut + Otiip)] - hu^ 



0. 



Noticing that ut = H ■ and = it • we may rewrite the above result as 



V{hu-k) - hil-x^n = 0. 



(3.89) 



Dividing by hu ■ k, we get 



V \n{-hu ■ k) - Vf] = 
u ■ k 



(3.90) 



Remark : The minus sign in the logarithm is justified by the fact, that since u and k are both future 
directed timelike vectors, u ■ k < 0. 



^For more general motions see § VI. E of [50]. 
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The integrability condition of Eq. (3.90) is eitlier (i) = const (rigid rotation) or (ii) tlie factor 
in front of Vfi is a function of fl: 



u ■ k 



F{n). 



(3.91) 



We sliall discuss tliese two cases below. But before proceeding, let us give an equivalent form of 
Eq. (3.90), based on the specific angular momentum 



(3.92) 



This quantity is (minus) the quotient of the two Bernoulli quantities discussed in § 3.4.2. It is 
thus conserved along any fluid line in a stationary and axisymmetric spacetime, irrespective of 
the fluid motion. Its name stems from the fact that I is the "angular momentum" (per unit mass) 
hu ■ X divided by the "energy" (per unit mass) —hu ■ ^. By combining (3.77) with (3.85) and 
= {n- uj)/N [Eq. (3.31)], we get 



At the Newtonian limit, B^l, N^l, uj^Q and I ~ Or^ sin^ 6. 

Since u ■ k = ut{l — 0,1), we may rewrite the equation of motion (3.90) as 

£ 



(3.93) 



vin[-hut{i-ne)] + 



1 - m 



vn = 0. 



(3.94) 



Besides, we have u^u^ = Utu^ + n^u''' = n*(uj + fin^) = u^Ut{l — n£), so that the normalization 



u^u^^ = — 1 leads to 



ut{i - ni) 



Therefore, Eq. (3.94) can be recast as 



Vln 



+ 



1 - n£ 



= 



(3.95) 



This is a variant of the equation of motion (3.90). 



Rigid rotation 

If = const, then Vr2 = and Eq. (3.90) leads immediately to the first integral of motion 

ln{—hu ■ k) = const. (3.96) 

Thanks to (3.86) and (3.31), we have 

N 

(3.97) 



u-k=[u) ^a,^— 

r/N 



so that the first integral (3.96) becomes 



H + u — InT = const . 



(3.98) 
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This first integral of motion has been obtained first in 1965 by Boyer [25] for incompressible fluids^ 
and generalized to compressible fluids by Hartle &: Sharp (1965) [59] and Boyer & Lindquist 
(1966) [27]. Note the difference of sign in front of InF with respect to the Bernoulli theorem 
(3.78). 

Remark : Another important difference is that (3.78) provides a constant along each fluid line, but 
this constant may vary from one fluid line to another one, whereas (3.98) involves a constant 
throughout the entire star. However, (3.98) requires the fluid to be in a pure circular motion, while 
the Bernoulli theorem (3.78) is valid for any stationary flow. 

The Newtonian limit of (3.98) is readily obtained since v ~ ^> and InT ~ C/V2 (cf. § 3.4.2); 
it reads 



if + ^> 



2 



const 



(3.99) 



Newt. 



Again notice the change of sign in front of f7^/2 with respect to the Bernoulli expression (3.83). 
Notice also that C/^ = ^^r^^\T?B [Eq. (3.41)] and that 



tot 



<& r^^r^sin^i 

2 



(3.100) 



is the total potential (i.e. that generating the gravitational force and the centrifugal one) in the 
frame rotating at the angular velocity ^2 about the rotation axis. 



Differential rotation 

Equation (3.90) gives, thanks to (3.79), (3.91) and (3.97), 

V (if + - InT) + F(J1)VJ7 = 0. 
We deduce immediately the first integral of motion : 



i7 + zy-lnr + / F{n')dn' = const 
Jo 



(3.101) 



Remark : The lower boundary in the above integral has been set to zero but it can be changed to any 
value: this will simply change the value of the constant in the right-hand side. 

Let us make explicit Eq. (3.91), by means of respectively Eqs. (3.97), (3.30), (3.31), (3.35) and 
(3.32): 



1 



N^[l - {B/N)^{n - u})^r^ sin2 ' 



Hence 



Fin) 



m - B'^{n-uj)'^r^ sin^i 



(3.102) 



^More precisely, the first integral obtained by Boyer is (within our notations) A'^^(l — U^) = const//*^. Taking 
the square root and the logarithm, we obtain (3.98). 
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At the Newtonian limit, tliis relation reduces to 



Newt. 



(3.103) 



and (3.101) to 



H + <^ + 

2 



/ F{n')dn' = const 
Jo 



(3.104) 



Newt. 



Equation (3.103) implies that must be a function of the distance from the rotation axis: 



Q = ^}{r sin ( 



Newt. 



(3.105) 



We thus recover a well known result for Newtonian differentially rotating stars, the so-called 
Poincare-Wavre theorem (see e.g. § 3.2.1 of Ref. [101]). In the general relativistic case, once 
some choice of the function F is made, one must solve Eq. (3.102) in Q to get the value of Q at 
each point. The first integral (3.101) has been first written^ by Boyer (1966) [2(3]. The relation 
(3.102) has been presented by Bonazzola & Machio (1971) [23] (see also the work by Abramowicz 
(1971, 1974) [1, 2]). 

From the alternative form (3.95) of the equation of fluid motion, we get, as an integrability 
condition, that — ili) must be a function of ft, i.e. that i itselft must be a function of Q: 



The relation between this function and F is obtained by comparing (3.94) with (3.90): 



F{n) 



1 - 



i + nF{n)' 



(3.106) 



(3.107) 



3.4.4 Stellar surface and maximum rotation velocity 

Let us consider a rigidly rotating star and rewrite the first integral of motion (3.98) as 

H = Hc + i^c-'^ + lrir, (3.108) 

where He and are the values of H and at r = [we have In Fc = since [7 = at r = by 
Eq. (3.32)]. The surface of the star is defined by the vanishing of the pressure: 



stellar surface: p = , 



(3.109) 



which expresses the local equilibrium with respect to the ambient vacuum. Let Hq be the value 
of H corresponding to p = (usually, Hq = 0). As a solution of the Poisson-like equation (3.14), 
u is an increasing function of r, from i^c < to = at r = +00. For instance, at the Newtonian 
limit and in spherical symmetry, v = —M/r outside the star. Therefore — is a decreasing 
function of r and if InF = (no rotation), (3.108) implies that H is a decreasing function of r as 
well. Provided that He > Hq, it always reaches Hq at some point, thereby defining the surface 
of the star. If the star is rotating, InF is an increasing function of r. It takes its maximum value 



''in the equivalent form N^(l - [/^) = [1 + r(n)f/h^, the link between the functions r(Sl) and F{Q,) being 

in[i + r(!^)] = - /" F(n')dn' 
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in the equatorial plane 9 = 7r/2. For instance, for a Newtonian star, InF = (f2^/2)r^ sin^ 0. If 
Q is too large, it could be that InF compensate the decay of —v in such a way that H remains 
always larger than Hq. Therefore, for a given value of H^, there exists a maximum rotation 
speed. Physically this corresponds to the fluid particles reaching the orbital Keplerian velocity 
at the stellar equator. For this reason, the maximal value of is called the Keplerian limit 
and denoted fix- 

3.5 Numerical resolution 

Having established all the relevant equations in the preceding sections, let us now discuss briefly 
their numerical resolution to get a model of a rotating star. 

3.5.1 The self-consistent-field method 

To be specific, we consider rigidly rotating stars: = const, but the following discussion can be 
adapted to take into account differential rotation. A standard algorithm to numerically construct 
a rotating stellar model is the following one. Choose 

• a barotropic EOS, of the form 

e = £{H) and p = p{H), (3.110) 
assuming that H = corresponds to p = 0; 

• some central value He of the log-enthalpy H; 

• some value for the constant angular velocity Q. 

In addition assume some values for the four metric functions N, lo, A and B. These may be 
crude values, such as the flat spacetime ones: N = A = B = 1 and a; = 0. Next, initialize U to 
zero and the enthalpy to a crude profile, like H = Hc{l — r"^ /B?) where R is some radius chosen, 
also crudely, around the expected radius of the final model. Then 

1. Via the EOS (3.110), evaluate e and p. The surface of the star is then defined by p = 0. 

2. From U , F, e and p, compute via Eqs. (3.36)-(3.38) the matter sources p^, S, 5*"^,, S^g, 
S'^ip that appear in the right-hand side of Einstein equations. 

3. Solve the Einstein equations (3.14)-(3.17) as flat Poisson-like equations, using the current 
value of v, co, A and B on their right-hand sides, thereby getting new values of = expi/, 
uj, A and B. 

4. Compute the velocity field from Eq. (3.32): U = (B /N){i} — uj)r sin 9; evaluate the Lorentz 
factor via Eq. (3.35): F = (1 - C/2)-V2. 

5. Use the first integral of fluid motion in the form (3.108): 

if = ifc + t'c - + In F 
to compute H in all space. Go to step 1. 
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In practice, this method leads to a unique solution for a given value of the input parameters 
{He, ri) and a fixed EOS. 

The above method is called the self- consistent- field method (SCF) and was introduced 
for computing Newtonian rotating stars by Ostriker &: Mark (1968) 1] (see Ref. [88] for some 
historical account). The first application of this method to compute rapidly rotating stellar 
models in full general relativity, along the lines sketched above is due to Bonazzola &; Maschio 
(1971) [23]. 

Remark : The procedure described above is only a sketch of the actual algorithm implemented in a 
numerical code; the latter is generally more complicated, involving some relaxation, progressive 
increase of from null value, convergence to a given mass, etc. 

The mathematical analysis of the self-consistent-field method applied to the Einstein equa- 
tions (3.14)-(3.17) has been performed by Schaudt & Pfister [92, 91, 86]. They proved the 
exponential convergence of the method and that the exterior and interior Dirichlet problem for 
relativistic rotating stars is solvable, under some conditions on the boundary data. The full 
mathematical demonstration of existence of rotating fluid stars in general relativity has been 
given by Heilig in 1995 [02], by a method different from the self-consistent-field, and only for 
sufficiently small angular velocities 0. Let us stress that, up to now, we do not know any exact 
solution describing both the interior and the exterior of a rotating star in general relativity (see 
e.g. [79, 93] for a discussion). 

3.5.2 Numerical codes 

The first numerical models of rotating relativistic stars have been computed by Hartle & Thorne 
(1968) [61] in the approximation of slow rotation (cf. § 3.2.5). The first rapidly rotating mod- 
els were obtained by Bonazzola &; Maschio (1971) [23] using a formulation very close to that 
presented here, namely based on Lewis-Papapetrou coordinates (the cylindrical version of QI 
coordinates, cf. § 2.2.2) and applying the self-consistent-field method to a system of equations 
similar to (3.14)-(3.17) (basically using the variable ^^—gu instead of N). All the numerical codes 
developed since then have been based on the self-consistent-field method. Almost all make use of 
QI coordinates, an exception being the code of Lin & Novak (2006) [77] based on Dirac gauge. In 
1972, Wilson [106] used the Bardeen- Wagoner formulation [14] (cf. § 3.2.5) to compute models 
of differentially rotating stars without an explicit EOS. In 1974 Bonazzola &: Schneider (1974) 
[24] improved Bonazzola &: Maschio's code and computed the first models models of rapidly 
rotating relativistic stars constructed on an explicit EOS (an ideal Fermi gas of neutron). In 
1975, Butterworth & Ipser [29, 30] developed a code based on the Bardeen- Wagoner formulation 
and computed homogeneous rotating bodies [29, 30] as well as polytropes [28]. Friedman, Ipser 
&; Parker (1986, 1989) [42, 43] applied their method to some realistic (i.e. based onto detailed 
microphysics calculations) EOS. This study has been extended by Lattimer et al. (1990) [70] to 
include more EOS. 

In 1989, Komatsu, Eriguchi &; Hachisu (KEH) [66] developed a new code, also based on the 
Bardeen- Wagoner formulation, and applied it to differentially rotating polytropes [67]. Their 
method has been improved by Cook, Shapiro & Teukolsky (GST) (1992,1994) [37, 38], by intro- 
ducing the radial variable s = r/{r -\- req), where r^q is the r coordinate of the stellar equator, 
in order that the computational domain extends to spatial infinity — the only place where the 
boundary conditions are known in advance [Eq. (3.23)]. In 1995, Stergioulas & Friedman [96] 
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Figure 3.1: Isocontours of the proper energy density e in tlie meridional plane = for the rotating 
neutron star model discussed in § 3.5.3. The coordinates {x, z) are defined by x := rsinO and z := rcosd. 
The thick solid line marks the surface of the star. 



developed their own code based on the CST scheme, improving the accuracy, leading to the 
famous public domain code rns [113]. 

In 1993, Bonazzola, Gourgoulhon, Salgado &: Marck (BGSM) [22] developed a new code based 
a formulation very close to that presented here (cf. § 3.2.5), using spectral methods [53], whereas 
all previous codes employed finite differences. This code has been improved in 1998 [21, 49] (by 
introducing numerical domains adapted to the surface of the star) and incorporated into the 
LORENE library, to become the public domain code nrotstar [112]. It is presented briefly in 
Appendix B. 

In 1998, Nozawa et al. [S3] have performed a detailed comparison of the KEH code, the rns 
code and the BGSM code, showing that the relative difference between rns and BGSM is of the 
order of 10~^ (at most 10~^ in some extreme cases). The relative difference between KEH and 
BGSM is only of the order 10~^ to 10~^. This is due to the finite computational domain used 
by KEH, implying approximate boundary conditions for some finite value of r. On the contrary, 
thanks to some compact ification, both rns and BGSM are based on a computational domain 
that extends to r = +oo where the exact boundary conditions (3.23) can be imposed. 

A new multi-domain spectral code has been developed by Ansorg, Kleinwaechter &: Meinel 
(AKM) in 2002 [8, 6]. It is extremely accurate, leading to rotating stellar models to machine 
accuracy (i. e. 10~^^ for 16-digits computations). It is has been used recently to compute 
sequences of differentially rotating stars, with a very high degree of differential rotating, leading 
to toroidal stellar shapes [7]. 

For a more extended discussion of the numerical codes, see the review article by Stergioulas 
[95]. In particular, a comparison between the AKM, KEH, BGSM, rns and Lorene/rotstar 
codes can be found in Table 2 of that article. 
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Gravitational mass M 


1.400 Mq 


Bar yon mass Mb 


1.542 M0 


Rotation rrequency iZ/(27rj 


716 Hz 


1 11 ill TT 

Central log-entnalpy He 


0.2262 c 


Central proper baryon density nt^c 


0.5301 fm~^ 


Central proper energy density 


5.7838 pnucc^ 


Central pressure Pc 


0.8628 pnucc^ 


Coordinate equatorial radius rgq 


9.867 km 


Coordinate polar radius 


8.649 km 


Axis ratio Tp/rcq 


0.8763 


Circumferential equat. radius i?circ 


12.08 km 


Compactness GM/(c /tcircj 


0.1711 


A 1 j_ T 

Angular momentum J 


0.7238 GMq/c 


Kerr parameter cJ/[GM ) 


0.3693 


Moment of inertia / 


1.417 X 10"** kgm 


Kinetic energy ratio T/W 


n HQ /I Q 


Velocity at the equator Ueq 


0.1967c 


Redshift from equator, forward z^^ 


-0.01362 


Redshift from equator, backward 


0.5529 


Redshift from pole Zp 


0.2618 



Table 3.1: Properties of the rotating neutron star model presented in § 3.5.3; 
Mq := 1.989 X 10^0 kg, p^uc := 1-66 x 10^^ kgm~^ 

3.5.3 An illustrative solution 

In order to provide some view of the various fields in a relativistic rotating star, we present 
here some plots from a numerical solution^ obtained by the Lorene/nrotstar code described in 
Appendix B. The solution corresponds to a neutron star of mass M = 1.4 Mq rigidly rotating 
at the frequency Q/{27r) = 716 Hz. This value has been chosen for it is the highest one among 
observed neutron stars: it is achieved by the pulsar PSR J1748-2446ad discovered in 2006 [63]. 
The EOS is the following one (see [•'),')] for details): 

• for the core: the model A18+(5u+UIX* of Akmal, Pandharipande & Ravenhall (1998) [3], 
describing a matter of neutrons, protons, electrons and muons via a Hamiltonian including 
two-body and three-body interactions, as well as relativistic corrections; 

• for the inner crust: the SLy4 model of Douchin & Haensel [41]; 

• for the outer crust: the Haensel &: Pichon (1994) model [54], which is based on the exper- 
imental masses of neutron rich nuclei. 

A global (coordinate-based) view of the star is provided in Fig. 3.1 and some of its global 
properties are listed in Table 3.1. All these quantities are properly defined in Chap. 4. The 
metric functions A^, oj and A are displayed in Figs. 3.2, 3.3 and 3.4 respectively. It appears that 

^This solution can be reproduced by using the parameter files stored in the directory 
Lorene/Codes/Nrotstar/Paraineters/GR/APR_l . 4Msol_716Hz 
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Figure 3.2: Profile of the lapse function N in two different directions: 6* = (rotation axis, dashed line) 
and 9 ^ tt/2 (equatorial plane, solid line). The vertical solid line (resp. dashed line) marks the location 
of the stellar surface in the equatorial plane (resp. along the rotation axis). 
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Figure 3.4: Same as Fig. 3.2 but for the metric coefficient A. 
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Figure 3.5: Difference between the metric coefficients B and A in the equatorial plane. The vertical 
solid line marks the location of the stellar surface in that plane. 
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6 




r [km] 

Figure 3.6: Profile of the fluid proper energy density e in two different directions: 6* = (rotation 
axis, dashed line) and 6 = 7r/2 (equatorial plane, solid line). The small vertical solid line (resp. dashed 
line) marks the location of the stellar surface in the equatorial plane (resp. along the rotation axis). 
Pnuc :=1.66x lO^kgni-^ 

the deviations from the flat space values (A^ = 1, a; = and A = 1) are maximal at the stellar 
center. One may also check on these figures the tendency to fulfill the boundary conditions (3.23) 
when r — )• +oo. Moreover, the faster decay of uj [in cf. Eq. (3.27)] with respect to that of 

A [in 1/r, cf. (3.63)] is clearly apparent. 

Remark : Figures 3.2 to 3.4 have been truncated at r = 30 km for graphical convenience, but the 
computational domain of the Lorene/nrotstar code extends to r = +oo. 

The metric function B turns out to be very close to A: on the rotation axis, it is exactly equal 
to A [cf. Eq. (2.36)] and, in the equatorial plane, the difference S — ^ is at most 2.3 x 10^'^, as 
shown in Fig. 3.5. 

Two fluid quantities are displayed in Figs. 3.6 and 3.7: the proper energy density e and the 
fluid velocity U with respect to the ZAMO [cf. Eqs. (3.32) and (3.33)]. In the first figure, note 
the presence of the neutron star crust, of width ~ 0.4 km for 6 = and ~ 0.6 km for 9 = 7r/2. 
The profile of U shown in Fig. 3.7 is deceptively simple: it looks very close to the Newtonian 
profile Or [Eq. (3.41) with 6 = vr/2], with one difference however: the mean slope is not Q but 
~ 1.3 0. This slope is approximately constant due to some numerical compensation of the two 
factors B/N and (1 —uj/Q) that appear in Eq. (3.32) once rewritten as U = {B /N)[l — uj /Vt) Vtr. 
Indeed B/N is a decreasing function of r, taking the values 2.36 at r — and 1.52 at r — r^q. 
On the contrary, 1 — w/O is an increasing function, taking the values 0.54 at r = and 0.87 at 
r = Tgq. This results in (B/N)(l — uj/Q) taking almost the same values r = and r = r^q (1.27 
and 1.32 respectively), hence the straight behavior of the plot of U{r). Notice that at the stellar 
surface, f/ is a non-negligible fraction of the velocity of light: U ~ 0.2 c. 
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Figure 3.7: Profile of the norm U of the fluid velocity with respect to the ZAMO, in the equatorial 
plane. 
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4.1 Total baryon number 

The baryon 4-current is 



3h = n^u, 



(4.1) 



where nb is the baryon number density in the fluid frame and w is the fluid 4- velocity. For a stellar 
model computed along the method sketched in § 3.5.1, nb is obtained from the log-enthalpy H 
and the EOS (3.110) via the formula 



nb 



e + p 



exp{-H), 



(4.2) 



which is easily deduced from Eqs. (3.79) and (3.73). The baryon 4-current obeys the conservation 
equation (3.70): 

= 0. (4.3) 
Thanks to this property, it can be shown, via the Gauss-Ostrogradski theorem, that the integral 



/ Jh ' dx^ dx^ dx^ 



(4. 
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is independent of the hypersurface S^, i.e. of the sHcing (St)jgiR of spacetime. It gives the total 
number of baryons in the star. Note that, by definition of the baryon 4-current, the integrand, 
= — jb • n, is nothing but the baryon number density as measured by the ZAMO, whose 4- 
velocity is n. Thanks to the relation F = — n • w [Eq. (3.34)], we have = Tn\^. Choosing 
(x^,x^,x^) = {r,9,ip), we have in addition ^77 = A^Br'^ sinO [Eq. (3.3)], so that expression (4.4) 
can be exphcited as 



Tn^A^Br^ sm 9 dr dO dip 



St 



(4.5) 



We may convert A/" to a mass, by multiplying it by the mean baryon mass mb — 1.66 X 10"^^ kg : 



Mb := mbAA 



(4.6) 



Mb is called the baryon mass of the star. 



4.2 Mass 

Due to the mass-energy equivalence, the mass of the star must involve all forms of energy, 
including the gravitational field energy. Now the latter is not localizable (as a consequence of 
the Equivalence Principle): there is no "energy density" of the gravitational field. This implies 
that the mass of the star cannot be defined as a integral over of some "mass density" as in 
(4.4), which gives the N as the integral of the baryon density = — jb ' n over S^. 

The are two major concepts of mass in general relativity (see Refs. [46], [100], [64] for more 
details) : the ADM mass which applies to any asymptotically flat spacetime and the Komar 
mass which holds only for stationary spacetimes. In the present case, both coincide. We shall 
examine first the Komar mass definition. 



4.2.1 Komar mass 

Given the Killing vector ^ associated with stationarity, the Komar mass is 



M 



d) V^e^ dS. 



■in 



(4.7) 



where ^ is any closed 2-surface (sphere) surrounding the star and dSfj^i, is the area-element 
2-form normal to J/". In the case where is entirely contained in the hypersurface St and is 
defined by r = const, it is spanned by the coordinates (9, ip) and 



dS^u = {s^n^ - s^n^)^ d9d(p, 



(4.8) 



where the n^'s are given by (1.54), the s^'s are the covariant components of the unit normal 
s of ^ in (St, 7) and q is the determinant of the components Qab of the metric induced by 7 
(or equivalently g) on From the diagonal form (2.39) of 7, we have s = s^dr and is 
determined by the normalization condition s-s = 1. Hence 



= (0, 0, 0) and s„ = (0, A, 0,0). 



(4.9) 



4.2 Mass 
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Since 5^ is defined by \t = const, r = const}, we read on (2.38) the following components of the 
metric q on S^: 

(4.10) 



Hence 



qah dx" dx^ = [A^dO^ + sin^ 9 dip^) . 
-y/g = ABr^ sin 6. 



(4.11) 



Remark : The Komar mass as defined by (4.7) is the flux integral through ^ of the type (2, 0) tensor 
V^. In dimension 3, there is locally only one direction s normal to a given sphere, so that instead 
of the normal 2-form (4.8), one considers the normal 1-form dSi ~ Siy^ dO dip; flux integrals are 
then defined for vectors, i.e. type (1,0) tensors. In the 4-dimensional spacetime (./#,g), there 
exists a whole plane of directions normal to the sphere and the vectors (n, s) which appear in (4.8) 
form a basis of that plane. 

A priori the quantity M defined by (4.7) should depend on the choice of the surface 
However, owing to the fact that ^ is a Killing vector, this is not the long as is located 

outside the star. This was shown by Komar in 1959 [()'>]■ To prove it, one may use the Gauss- 
Ostrogradski theorem, the Killing equation (2.4) and the Einstein equation (1.24) to rewrite the 
integral (4.7) as"*^ 

M = 2 / T{fi,i) -]-Tn- i ^drdOdip, (4.12) 

where T is the energy-momentum tensor and T := g^'^T^^j its trace with respect to g. From 
(4.12), it is clear that outside the support of T, i.e. outside the star, there is no contribution to 
M, hence the independency of M upon the choice of 5^. 

Expressing (4.7) in terms of 3+1 quantities, we get (cf. § 7.6.2 of [±u] for details). 



47r 



.9" 



1.13) 



In the present case, substituting (4.9) for s\ (2.43)-(2.44) for Kij, (2.37) for /3-'' and (4.11) for 
y^, we have 



M 



1 

47r 



.9' 



dN B'^r'^ sin2 6 du 
'W 2N ^'d^ 



Br^ sin 6 d9 dip. 



(4.14) 



If we let the radius of the sphere extend to +cxd, then the term containing u does no longer 
contribute to the integral, since (3.27) implies 



dr 

Taking into account i? — t- 1 when r — t- +oo, we may then write 



1 



M 



— lim 

47r y^oo 



dN 2 

— — r svnOdOdip 



(4.15) 



where the limit means that the integral is to be taken on a sphere =y of radius r — )• +oo. Plugging 
the asymptotic behavior (3.26) of N into the above formula yields Kq = —M, hence 



N = 1 



(4.16) 



^See e.g. § 7.6.1 of [4()] for the detailed demonstration. 
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Let us now turn to the volume expression of M provided by formula (4.12). Introducing in it 
the 3+1 decomposition (1.50) of the energy-momentum tensor, and making use of ^ = Nix + (3 
[Eq. (1.51)] with dt = ^, we get 



M=j [N{E + S)-2pip']^drded^. 



1.17) 



In the present case, taking into account expression (3.37) for pi, (2.37) for /3* and (3.3) for ^77, 
we obtain 



M = [ [N{E + S) + 2ujB{E + p)Ursme]A'^Br'^smedrd9dcp 



At the Newtonian hmit, only the term in ~ /) remains in this expression: 



M = j pr sin 6 dr d9 dip 



(4.18) 



(4.19) 



Newt. 



Of course, we recognize the expression of the mass as the integral of the mass density p. 



4.2.2 ADM mass 

The ADM mass has been introduced by Arnowitt, Deser & Misner in 1962 [9] for any asymp- 
totically flat spacelike hypersurface S^. It is expressible as [46, 64] 



Madm = 7^ lim 

Ibvr .y-+oo jy 



P^7i, -A(/%0 s'^dOd^ 



where (i) / is a fiducial flat 3-metric on Sj: fij = diag(l, r , r sin 6) and /*■ 
and (ii) T) stands for the connection associated with /. 



diag(l, r 



(4.20) 



-2 —2 ■ - 

,r sm 



Remark : Contrary to the Komar mass definition (4.7), the surface integral in (4.20) has to be taken 
at infinity. Note also that the ADM mass does not require stationarity to be defined. In non- 
stationary spacetimes, the Komar mass can still be defined by an integral of the form (4.7) where 
^ is an asymptotic Killing vector, but then the surface has to be located at infinity. 

We leave as an exercise to the reader to show that with the form (2.39) of 7jj, the formula 
(4.20) can be explicited as 



Ibvr y-^oo jy 



— {A^ + B^) + 

or r 



r sin dO dip 



(4.21) 



Although the two definitions look very different, Beig (1978) [16] and Ashtekar & Magnon- 
Ashtekar (1979) [10] have shown that for a stationary spacetime, the ADM mass is always equal 
to the Komar mass, provided that Et is orthogonal to the Killing vector ^ at spatial infinity, 
which is our case here: 

(4.22) 



Madm = M 



Therefore, for a stationary rotating star, the Komar mass and ADM mass coincide. This 
common value, denoted here by M, is sometimes called gravitational mass, to distinguish it 
from the baryon mass defined by (4.6). 
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4.2.3 Binding energy 

The binding energy of the star is defined by as the difference between the gravitational mass 
and baryon mass: 



^bind := Af - Mb . (4.23) 



Physically, l-Ebindl is the energy one should spend to disperse all the particles (baryons) consti- 
tuting the star to infinity. The star is bound iff E'bind < 0. If -Ebind ^ 0, the star would be 
unstable and would explode, releasing the energy -Ebind- 



4.3 Virial identities 
4.3.1 GRV3 

Taking the identity of ADM and Komar mass as a starting point, Gourgoulhon &: Bonazzola 
(1994) [48] have derived a formula that generalizes the Newtonian virial theorem to general 
relativity for stationary spacetimes. It is called GRV3, for General Relativistic Virial theorem^ 
the "3" recalling that the integral is a 3-dimensional one. For a rotating star, it reads [Eq. (43) 
of [48]] 



AttS - 4t (dvdv \-dAdB] + — ( 4t - 4j 

\ 2AB J 2r\A^ 



1 fdA 1 dA 
A\dr r tan 6 89 



1 fdB 1 dB 

+ 



2B \dr r tan 9 89 
with S given by (3.38) : 



o d2 oj-qS q 

+ g^2^2jv2 ^^^^ \A^Br'^s\n9dTd9dip = 0, (4.24) 



S = 2,p + {E + p)U'^ . 

At the Newtonian limit, S ~ 3p + pU"^, v ~ ^>/c^. In A ~ -<^/c^, Ini? ~ -^jc? , A ~ 1, ~ 1 
and a; = 0, so that, after division by 47r, (4.24) becomes 



2^kin + 3P + ^grav = (4.25) 
I Newt. 



where 

EY,n = \\ pU^ r'^ sin 9 drd9 dip (4.26) 

P= pr^sm9drd9dip (4.27) 

-Egrav = ^/ {V<!>f r'^ sine drd9 dip. (4.28) 

E'kin and -Egrav ai'e respectively the kinetic energy and the gravitational potential energy of the 
star, and we recognize in (4.25) the classical virial theorem for stationary configurations (see e.g. 
Eq. (2.180) of Ref. [101]). 

The GRV3 identity proved to be very useful to check numerical solutions [83, 95]. 
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4.3.2 GRV2 



Another virial identity has been obtained by Bonazzola (1973) [19] (see also Ref. [22, 20]). 
Following [20], it is called GRV2, because contrary to GRV3, it involves a 2-dimensional integral. 
The derivation of GRV2 is quite straightforward once we consider the last equation of the Einstein 
system (3.14)-(3.17). The operator on the left-hand side of (3.17) is the 2-dimensional flat 
Laplacian A2. It operates in the plane spanned by the polar coordinates {r,6) and its Green 
function is 

G{x,x') = ^ln\\x-x'\\, (4.29) 
where x denotes the vector of of coordinates (x") = {r,6) and \\x — x'\\ is the norm of x — x' 



with respect to the flat Euclidean metric of 
assuming its right-hand side is given, is 



Accordingly, the generic solution of (3.17), 



In A + 1^ = — Inr + O ( - 
2tt \r 



where / is the integral of the right-hand side of (3.17), 



a := 87r^2 w + 



352r2 sin2 e 
4/V2 



dudco — dudu. 



over with respect to the flat area element r dr dO: 



I 



r=+oo /•9=2iT 



r=0 



ar dr dO. 



(4.30) 



(4.31) 



(4.32) 



In the above integral, the coordinate 9 spans the interval [0, 27r], whereas the QI coordinate 
runs in [0, tt] only. This is because the meridional surface Mtip covers only half of M^. To define 
(4.32), we therefore use the analytical continuation of a from [0, tt] to [vr, 27r]. More precisely, since 
the regularity conditions imply that a can be expanded into a series of cos{k6), the analytical 
continuation satisfies 

V0G[^,2^], a{r,9) = a{r,2Tr-9). (4.33) 

Now, when r — )• -|-oo, the asymptotic behavior of ly is 0(1 /r) [cf. Eq. (3.26)] and the same things 
holds for In A [see (3.63) for the spherically symmetric case]. To have a regular asymptotically 
flat solution, we must therefore have / = in (4.30). In view of (4.31), (4.33) and expression 
(3.38) for iS"^<^, this implies 



r=0 



9=0 



8TrA'^\p + iE + p)U'^] + 



3B^r^ sin^ 
4iV2 



dujduj — dvdv > rdrd6 = 



(4.34) 



This is the GRV2 identity. Note that the integration domain in has been reduced to [0, vr] 
thanks to the property (4.33), thanks to which / is twice the integral in the left-hand side of 
(4.34). As GRV3, it is very useful to check the accuracy of numerical solutions [83, 95]. 
The Newtonian limit of GRV2 is 



/•r=+oo /•9=iT 






Jr=o Ve=o 


r dr d9 = 



(4.35) 



Newt. 
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The identity (4.35) has been exhibited in 1939 by Chandrasekhar [35] in the special case of 
spherical symmetry, but it seems that it has not been considered for axisymmetric Newtonian 
bodies before the work [20]. 

Remark : Although it also involves a 2-dimensional Laplaciaii, Eq. (3.16) of the Einstein system does 
not give birth to a virial identity because of the sin 6 factor on its right-hand side, which leads to 
a{r, 6) = —(j(r, 27r — 9) instead of (4.33), implying that / is automatically zero. 



4.4 Angular momentum 

The Komar angular momentum (or angular momentum for brevety) is defined by a 
integral similar to that used to define the Komar mass in § 4.2.1 but replacing the stationarity 
Killing vector ^ by the axisymmetry Killing vector X'- 



1.36) 



where ^ is any closed 2-surface (sphere) surrounding the star and dS^v is the area-element 
2-form normal to J/" given by (4.8). As for the Komar mass, J is independent of the choice of 
^ provided the star is entirely enclosed in =5^. 

In terms of the 3+1 quantities J is expressible as (cf. [Ki] for details) 



^in ^ij^'^^V^de dtp. 



(4.37) 



By means of the Gauss-Ostrogradski theorem and Einstein equation, we can get a volume ex- 
pression of J similar to (4.12), differing only by the replacement of ^ by % and the overall factor 
2 by —1. Since x is tangent to S^, we have T{n,x) = —p- X ^ ■ X = 0, so that the final 
expression reduces to 

J= pix'^drdOdif. (4.38) 
For a perfect fiuid rotating star, piX^ = p^p and we get, by combining Eqs. (3.37) and (3.3), 



J = I {E + p)UA^B'^r'^sm^edrdedip 



At the Newtonian limit, this expression reduces to 



J = I pUr sin 6 x r sin 9 dr dO dip 

St 



Newt. 



(4.39) 



(4.40) 



We recognize the z-component of the star total angular momentum. 

Let us make explicit the surface integral expression of J, Eq. (4.37), in the present case. We 
have from Eq. (2.43)-(2.44) and (4.9), 



B'^r'^ SVC? 6 doj o /i^?"^ i 

— ~ — r smt^-— when r — )• +oo. 

dr dr 
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Therefore taking the hmit ^ — t- oo in Eq. (4.37) gives 

J - 



lim </ '^^r^ siv? 9 d9 dip . 

\Qtt y^oo jy or 

Now, from the asymptotic behavior of lo as given by (3.27), du/dr = — 3-fCi/r^ + 0(l/r^). Hence 



J 



3K 



sin-^ 6 de 



3Ki 



10 o JO 

We may therefore rewrite the asymptotic behavior of uj as 



/ {i-cos^e)smede = 

Jo 2 



UJ 



'4 + 0^' 



(4,41) 



In particular, we recover the fact that if the star is not rotating (J = 0), then u = [cf. (2.42)]. 
Having defined both J and M, we can form the dimensionless ratio 



J 

IP 



(4.42) 



j is called the Kerr parameter. For a Kerr black hole of parameters (a, M), j = a/M and 
varies between (Schwarzschild black hole) and 1 (extreme Kerr black hole). 
Next, for a rigidly rotating star, one defines the moment of inertia by 



J 



(4.43) 



Remark : Due to the so-called supertranslation ambiguity, there is no equivalent of the ADM mass for 
the angular momentum (see e.g. [(i4] for a discussion). 



4.5 Circumferential radius 

The stellar equator is the closed line at the surface of the star [p = 0, cf. (3.109)] defined by 
t = const and 9 = tt/2 (equatorial plane). It has a constant value of the coordinate r, that we 
have already denoted by r^q. A coordinate-independent characterization of the stellar equator is 
the circumferential radius: 



Rc 



2^ 



(4.44) 



where ^ is the circumference of the star in the equatorial plane, i.e. the length of the equator 
as given by the metric tensor. According to the line element (2.38), we have 



1 

2^ 



r = rcq 
= 7r/2 



1 

2^ 



= 7r/2 



27r 



hence 



Rc 



5(req,7r/2)r, 



eq 



Of course, at the Newtonian limit, B = \ and R^ 



' eq- 



(4.45) 

For a relativistic star, we have 



B > 1 (cf. Fig. 3.4, taking into account that B is very close to A, as shown in Fig. 3.5), so that 

Reive ^ ^e 



' eq- 
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4.6.1 General expression of the redshift 

Consider a photon emitted by the fluid at the surface of the star. The photon follows a null 
geodesic ^ of spacetime and is described by its 4-momentum^ p, which is a future directed null 
vector: 

p-p = 0. (4.46) 
p is tangent to ^ in each point and obeys the geodesic equation 

Vp^p = 0. (4.47) 

For an observer who is receiving the photon, the relative shift in wavelength, called redshift, 

is 



A 



z :- 



rec ^>em 



Ac 



A. 



where Arcc is the photon's wavelength as measured by the observer (the "receiver") and Acm is the 
photon's wavelength in the fluid frame (the "emitter"). In terms of the photon energy E = hc/X, 
the above formula can be written 



Ere 



1. 



(4.49) 



Now, by the very definition of the photon 4-momentum, the photon energy measured by an 
observer O is 

Eo = -uo ■ p, (4.50) 

where uo is the 4- velocity of O. In the present case, Ucm = u (the fluid 4- velocity), whereas, for 
an observer at infinity and at rest with respect to the star u^cc = ^- Hence 



(4.51) 




where the indices 'em' and 'rec' indicate the location where the scalar products u ■ p and ^ • p 
are to be taken. 

A fundamental property for the scalar product ^ • p is to be conserved along the photon 
geodesic worldline. To show it, let us take the scalar product of (4.47) by ^; we have successively 

e.p^v^p^ = 0, 

p''v^(e.p")-pVv^e. = o, 



P^V,.{^uPn = 0, (4.52) 



■^We are using the same letter p as for the matter momentum density with respect to the ZAMO inside the 
star [cf. Eq. (1.32)], but no confusion should arise. 
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where we have used the Kilhng equation (2.4) to set to zero the term p^p'^V^^^: 



^ ^ ^ — — ' 



0. 



Equation (4.52) can be rewritten as 



VpK^ • p) = 



Since p is tangent to the photon worldhne, this relation shows that 



^ ■ p = const along the photon worldline 



(4.53) 
(4.54) 



Hence the symmetry encoded by ^ (stationarity) gives birth to a conserved quantity. This is 
another manifestation of the Noether theorem (cf. § 3.4.2). Similarly, we have 



X ' p = const along the photon worldline 



(4.55) 



The conservation laws (4.54) and (4.55) are particularly useful when we specialize u in (4.51) 
to the circular form (2.25): u = {T/N){^ + 17%) (we have used (3.31) to replace u* by T/N). We 
get then 

u-p= — ($^-p + nx-pj , 



N 



so that Eq. (4.51) becomes 



Thanks to (4.54), $, ■ p 



N 





p 


cm 




p 


rcc 



^■P 



(4.56) 



hence 



r 

N 



1 + 

cm \ ^ ' P 



(4.57) 



Note that we have suppressed the indices 'em' or 'rec' on x ' P ^-nd ^ • p thanks to (4.54) and 
(4.55). 

Remark : For difTerentially rotating stars, ft is not a constant and the value to be used in Eq. (4.57) is 
the value at the emission point. 



4.6.2 Polar redshift 

Let us apply formula (4.57) to an emission taking place at one of the poles of the star, i.e. at 
one of the two points of intersection of the stellar surface with the rotation axis: ^ = or ^ = tt. 
On the rotation axis X = [Eq. (2.6)]. Moreover, F = 1 since {/ = on the rotation axis [cf 
Eq. (3.32) with sinO = 0]. Accordingly, formula (4.57) simplifies to 



— - 1 

Np 



(4.58) 



where the index 'p' stands for polar values. 
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4.6.3 Equatorial redshifts 



Let us now consider the emission from the equator and in a direction tangent to it, i.e. an 
emission from the left or right edge of the star as seen by the distant observer. The photon 
4-momentum at the emission is then 



(4.59) 



with p'^ > (resp. p"^ < 0) for an emission in the forward (resp. backward) direction with 
respect to the fluid motion. The term in factor of Q in formula (4.57) is then 



X-P _X-^ + ax-X _ gt^ + ag^^ 



with 



a 



(4.60) 



^ • P i-i + Oii-X 9tt + agt^ p' 

Note that since p is a future directed null vector, we have > 0; in particular ^ 0, so that a 
is always well defined. 

Besides, Eq. (4.46) writes gtt{p*f + '^gt^P^P'^ + QwiP^^f = 0' i-^- 

g^^a^ + 2gt^a + gu = 0. (4.61) 
From this relation, we can write (4.60) as 

X P ^ 

Finally, we compute a by solving the second-order equation (4.61). We get two solutions: 



(4.62) 



9t^ 



gttgw Bruj ± N 



gifiifi Br 

where the second equality follows from the expression (2.34) of the metric components, specialized 
at 9 = Tr/2 (equatorial plane). Substituting the above value for a into Eq. (4.62) leads to 

Br 



XP 



i p - Bruj' 



(4.63) 



Since ^ and p are both future directed, we have always $ • p < 0. For a emission in the forward 
direction X P = B'^r'^ijp'^ — up'') > 0, so the — sign must selected in the above expression. On 
the contrary, the + sign corresponds to an emission in the backward direction. Inserting this 
result into Eq. (4.57), we thus get, for an emission in the forward direction 



Br 



N + Bruj 



1 



N + Broj 



1 



u 



OJ 



1 



1 



1 - U 



N + Bruj Vl - C/2 



1, 



I.e. 



1 



'^'^ N + Bru \ l + U 



1 - U 



(4.64) 



In the right-hand side, all the functions are to be taken at the emission point. Similarly choosing 
the -|- sign in (4.63) leads to the redshift in the backward direction: 



1 



N - Bru V 1-U 



l + U 



(4.65) 
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Remark : As a check of the above expressions, we may consider the flat spacetime hmit: N = 1 and 
w = 0. We then recover the standard Doppler formulas 

" V 1 + c/ 

They imply < 0, i.e. a blueshift (the emitter is moving towards the observer in the forward case) 
and > 0, i.e. a true redshift (the emitter is moving away from the observer in the backward 
case). 



1 + U 
1 - U 



4.7 Orbits around the star 
4.7.1 General properties of orbits 

Let us consider a test particle ^ of mass m <^ M orbiting the star. Let v be the 4-velocity of 
its components with respect to QI coordinates are 

^ f dt dr dO d(p\ , , 

where r is the proper time of The 4-momentum of the particle is 

p = mv. (4-67) 

The 4-velocity normalization relation -j; • -j; = — 1 is equivalent to 

pp = —m^. (4.68) 

Assuming that ^ is submitted only to the gravitation of the star, its worldline is a timelike 
geodesic of (^,gi). In particular = or equivalently 

Vjyp = 0. (4.69) 

Due to the symmetries of (^,gi), we have the same conserved quantities as those given by 
Eqs. (4.54)-(4.55) in the photon case: 



E := -i ■ p = -pt = const (4.70) 



L ■= X ' P = Pip = const . (4.71) 



Indeed the demonstration of (4.54)-(4.55) [cf. Eq. (4.52)] used only (i) the Killing character of ^ 
and X ^iid (ii) the property (4.69), which holds for a photon as well as for a massive particle. If 
the particle reaches spacelike infinity, where ^ n (the ZAMO 4-velocity), the constant E is the 
particle energy as measured by the asymptotic ZAMO Similarly, L is the angular momentum of 
the particle with respect to the rotation axis as measured by the same observer: when r — >■ -|-oo, 
L ~ mr'^ sin^ 6 dip/dt (assuming that asymptotically, the proper time of the ^ coincides with t, 
i.e. that ^ is not relativistic). 

Everywhere is space, the energy of ^ as measured by the ZAMO is [cf. Eq. (2.41)] 

^ZAMO = -n p = (^i+(^Xj "P = {^-i-p-^X-P 
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hence 

(4.72) 



Remark : Since 1 and uj ^ when r — > +oo, it is clear on this formula that E = Ezamo h the 

particle reaches spatial infinity. 



4.7.2 Orbits in the equatorial plane 



We restrict ourselves to orbits in the equatorial plane: = tt/2. Then, according to Eq. (4.66) 
and (4.67), = mdO/dr = 0. Given the form (2.32) of the metric coefficients, this yields 



0. 



Combining with (4.70)-(4.7f ), we conclude that we have three constants among the four compo- 
nents Pa- pt, Pe and p^. The fourth component, p,., is obtained from the normalization relation 
(4.68): 



g^'^PiiPu 



Substituting the values of g^^ given by (2.35), we get 



Pt _ r, ^PtPv> ,P^, Pe 
'iV2 iV2 ^ r-^A^ 



+ 



1 



Since 6 = ir/l, pt = —E, pe = and = L, we get : 

Pr ^ f IP ^^2 2 ^ 

A^ B^r^ 



Now, from (2.34), (4.66) and (4.67), 



g^^pf^ = A^p'' = A^mv'^ = A^m'^. 



Thus, we can rewrite (4.73) as 



1(1) +V(.fi.l) = 



where E and L are the constants of motion defined by 

- E - 
E := — and L :- 

m 

and 



L 



V{r,E,L) :-- 



1 



2.42 



(4.73) 



(4.74) 



(4.75) 



1.76) 



In view of Eq. (4.74), we conclude that the geodesic motion of a test particle in the equatorial 
plane is equivalent to the Newtonian motion of a particle in a one-dimensional space (spanned 
by the coordinate r) under the effective potential V{r,E,L). 
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4.7.3 Circular orbits 

Circular orbits are defined by r = const, i.e. by 



dr 
(fr 



0. 



(4.77) 
(4.78) 



Remark : Condition (4.77) is not sufficient to ensure a circular orbit, for it is satisfied at tlie periastron 
or apoastron of non-circular orbits. 

From Eq. (4.74), the two conditions (4.77)-(4.78) are equivalent to V = and dV /dr = 0, i.e. to 

-2 



E-u}L 







iV2 



, ^ v^dv -doj 
{E - ojL)— + L— 
or or 



L2 /d(3 1 



]^2^2 \ Qr r 



where 



i':=\nN and /? := InB. 



(4.79) 
(4.80) 

(4.81) 



As we did for the fluid, let us introduce the particle angular velocity as seen by a distant 
observer, 



dip 
dt 



(4.82) 



as well as the particle velocity V measured by the ZAMO. Then, similarly to Eqs. (2.25), (3.30) 
and (3.30) we have 

v = v' + n^x) , (4.83) 

(4.84) 

^ ^ ^ (4.85) 



V = T^ [fi + V] , 

T^ = Nv* and 



As for (3.33), the norm of V is \V\ with [cf. Eq. (3.32) with sin 6' = 1] 

B 



V := j^i^^ -u)r. 
The normalization relation i> • -i; = — 1 is then equivalent to 

From Eq. (4.70) and (4.83), we have 

E = -vt = -gtf.v'' = (iV^ - B^uj'^r'^y + ioB'^r'^v^n,^ = [N^ + B'^ur'^{n^ - to)] . 
Thanks to (4.85) and (4.86), we can write 



(4.86) 



(4.87) 



E = r,^{N + BujrV) 



(4i 
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Next, from Eq. (4.71) and (4.83), we have 

L = v^= g^^v^ = -uB'^r'^v^ + B^r'^v^n^ = v^B^r'^{n^ - w), 
i.e. thanks to (4.85) and (4.86), 



L = BrTcmV 



Combining (4.88) and (4.89), we have 



E-ujL = T^N 



(4.89) 



(4.90) 



This result leads to the following expression for the energy of as measured by the ZAMO and 
given by Eq. (4.72): 

' ' (4.91) 



-EzAMO = r,^ m 



We thus recover Einstein's mass-energy relation. This is not surprising since is the Lorentz 
factor of ^ with respect to the ZAMO. 

In view of the above relations, let us go back to the system (4.79)-(4.80). Thanks to (4.90) 
and (4.89), Eq. (4.79) is equivalent to 1 - r|, + T%,V'^ = 0. It is trivially satisfied since and 
V are related by (4.87). Inserting (4.90) and (4.89) into Eq. (4.80) gives rise to a second order 
equation for V: 



dp 



Br du 



dv 



^ + _M/^ V = 0. 

dr r J N dr dr 

This equation was first written by Bardeen (1970) [II]. It admits two solutions: 



(4.92) 



V± 



ur QUI _|_ 
N dr ^ 



TP' \dr) ^dr \ dr ^ r 



2 ( M . 1 

dr ^ r 



(4.93) 



This formula gives the circular orbit velocity with respect to the ZAMO as a function of r. All 
the metric coefficients in the right-hand side are to be taken at ^ = 7r/2. The -|- (resp. — ) sign 
corresponds to a motion in the same sens (resp. opposite sens) than the stellar rotation: the 
orbit is then called direct (resp. retrograde). 

At the Newtonian limit, a; — )■ 0, \df3/dr\ <ti 1/r and u ^/c^, so that (4.93) reduces to 



V± 



c'r or 



I.e. 



(4.94) 



Newt. 



We thus recover the standard formula obtained by equating the gravitational acceleration, — 



and the centrifugal one, —V'^/r. In particular, in spherical symmetry, <I> 
yields to the well known result : V± = ±y^GM/r. 



-GM/r and (4.94) 



Remark : At the Newtonian limit, = — V+, whereas in the relativistic case, if the star is rotating 

^ 0), y_ ^ ~v+. 
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4.7.4 Innermost stable circular orbit (ISCO) 



We have seen from Eq. (4.74) that a circular orbit corresponds to an extremum of the effective 
potential: dV/dr = 0. It is stable if, and only if, this extremum is actually a minimum, i.e. if, 
and only if, 

^ > 0. (4.95) 

Taking into account that V = and dV/dr = on a circular orbit, the above criteria is equivalent 
to 



E-uL 



N2 



+ 



{E - coL) 



Qj.2 



\ dr 



9^/3 4 dp 
_| — _|_ 2 

Qj.2 J. Qy. 



dr 



+ L 



+ 



Qy.2 



dr dr , 



iV2 V ar ' ^ 



(4.96) 



Substituting (4.90) for E — ujL and (4.89) for L, we get, after division by 



\ dr 



+ 



VBr 

IT 



d^u 

Qy.2 



Qj.2 r dr \ dr 



^du dio 
dr dr 



iV2 



§^1 >o. 

dr 



(4.97) 



The marginally stable orbit, if any, corresponds to the vanishing of the left-hand side in the above 
inequality. This orbit is called the innermost stable circular orbit (ISCO). To compute its 
location, one has to solve the equation formed by replacing the > sign in (4.97) by an equal sign: 



8ih 
dr' 



V dr 



+ 



VBr 



N \ dr"^ 



Qj.2 r dr \ dr 



^dv dio 
dr dr 



iV2 



duo 
dr 



^ =0. 



(4.98) 



In this equation, V is considered as the function of the metric potentials A^, B and oj given by 
(4.93) and all the potentials are to be taken at 9 = tt/2. So, once the metric is known, Eq. (4.98) 
is an equation in r. If it admits a solution r > rcq (i.e. outside the star), then the ISCO exists. 
Otherwise, all the circular orbits are stable down to the stellar surface. 

At the Newtonian limit (zv ~ <I>/c^, |/3| <C 1 and oj — ^ 0), Eq. (4.98) reduces to 



52^> 3 d^ 

Q^2 y. Qj. 



0, 



(4.99) 



where we have used (4.94) to express V"^. 

There are basically two causes for the existence of the ISCO: 

• The star is highly relativistic: for instance in spherical symmetry the metric outside the 
star is Schwarzschild metric (cf. § 3.3.3); if the star is compact enough, its surface is located 
under the famous f = 6M which defines the ISCO in Schwarzschild metric. Note that via 
the relation (3.59) between the areal coordinate f and the isotropic coordinate r, f = 6M 
corresponds to r = (5/2 + \/6)Af • 
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• The star deviates significantly from spherical symmetry: in Newtonian gravity, all the orbits 
are stable around a spherically symmetric star. But if the gravitational potential deviates 
significantly from $ = —GM/r, an ISCO may exist, even in Newtonian gravity. This was 
shown explicitly for very rapidly rotating low-mass strange quark stars in Ref. [109]. For 
instance, assuming that the deviation from spherical symmetry arises only from the mass 
quadrupole moment of the star Q (this is the leading non-spherical term when r — ?• -|-oo) , 



where P2{x) = (3x^ — l)/2 is the Legendre polynomial of degree 2. Inserting the above 
expression of <I> in Eq. (4.99), setting 9 = Tr/2 and solving for r leads to 




^P,(cos^) 



(4.100) 




(4.101) 



^Equation (4.100) is just Eq. (3.25) truncated at ^ = 2 and with equatorial symmetry, so that the term £ = 1 
vanishes. 
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A.l Lie derivative of a vector field 
A. 1.1 Introduction 

We have seen in § 1.1.5 that the definition of the derivative of a vector field v on ^ requires 
some extra structure on the manifold This can be some connection V (such as the Levi- 
Civita connection associated with the metric tensor g), leading to the covariant derivative Wv. 
Another extra structure can be some reference vector field iZ, leading to the concept of derivative 
of V along u] this is the Lie derivative presented in this Appendix. These two types of derivative 
generalize straightforwardly to any kind of tensor field^. 

A. 1.2 Definition 

Consider a vector field u on called hereafter the flow. Let v be another vector field on 
the variation of which is to be studied. We can use the flow u to transport the vector v from one 
point p to a neighbouring one q and then define rigorously the variation of v as the difference 
between the actual value of -y at g and the transported value via u. More precisely the definition 
of the Lie derivative of v with respect to u is as follows (see Fig. A.l). We first define the image 
$e(p) of the point p by the transport by an infinitesimal "distance" e along the field lines of u 
as ^e{p) = 9) where q is the point close to p such that p^ = eu{p). Besides, if we multiply the 
vector v{p) by some infinitesimal parameter A, it becomes an infinitesimal vector at p. Then 
there exists a unique point p' close to p such that Xv{p) = pf)'. We may transport the point p' 

^For the specific Itind of tensor fields constituted by antisymmetric multilinear forms, there exists a third type 
of derivative, which does not require any extra structure on : the exterior derivative (see the classical textbooks 
[82, 104, 98] or Ref. [45] for an introduction). 
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y(p) 



Figure A.l: Geometrical construction of tlie Lie derivative of a vector field: given a small pa- 
rameter A, each extremity of the arrow Xv is dragged by some small parameter e along u, to 
form the vector denoted by <^e(A'y). The latter is then compared with the actual value of Xv at 
the point q, the difference (divided by Ae) defining the Lie derivative £^ v. 



to a point q' along the field lines of it by the same "distance" e as that used to transport p to q: 
q' = ^eip') (see Fig. A.l). q^' is then an infinitesimal vector at q and we define the transport by 
the distance e of the vector v{p) along the field lines of u according to 



1 



X 



(A.l) 



^i;{v{p)) is vector in J3^(^). We may then subtract it from the actual value of the field v at q 
and define the Lie derivative of v along u by 



C^v:= lim - [v{q) - ^e{v{p))] 
e 



(A.2) 



If we consider a coordinate system (x") adapted to the field u in the sense that u = Oq where 
Sq is the first vector of the natural basis associated with the coordinates (x") (cf. § 1.1.1), then 
the Lie derivative is simply given by the partial derivative of the vector components with respect 



to a;^ 



dv°' 



In an arbitrary coordinate system, this formula is generalized to 



(A.3) 



(A.4) 



where use has been made of the standard notation >C^w" := (C^v)". The above relation shows 
that the Lie derivative of a vector with respect to another one is nothing but the commutator of 
these two vectors: 

■ v = \u. v] I. fA.5) 
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A. 2 Generalization to any tensor field 

The Lie derivative is extended to any tensor field by (i) demanding that for a scalar field /, 
f = df[u) and (ii) using the Leibniz rule. As a result, the Lie derivative T of a tensor 

field T of type is a tensor field of the same type, the components of which with respect to 

a given coordinate system (x") are 

^ _ rpai...ak _sr^ ai...a...ak du"' ST^ rpai...ak ^^"^ 

^u-^ /3i.../3f - " g^i,^ /3i.../3f /3i..A Q^a ^ /3i....t.../3^ g^/J^ 

i=l i=l 



In particular, for a 1-form, 

^x^^ ^ dx 

and, for a bilinear form, 



(A.6) 



^u^a = U^'^+UJ^— (A.7) 



^.^^.. = -^^+^..^ + ^«.^- (A.8) 

Notice that, thanks to the symmetry properties of the Christoffel symbols defined by (1.19), the 
partial derivatives the above equations can be replaced by the covariant derivatives V associated 
with the metric gr, yielding 

i=l 1=1 \ 

(A.9) 

Cu UJa = uf'V^UJa + ^iNo.U^', (A. 10) 

CuTap = u^'^^T^p + T^pVaU^" + T^^Npu^'. (A.ll) 
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Appendix B 

Lorene/nrotstar code 



The nrotstar code is a free software based on the C++ LORENE hbrary [112]. It comes along 
when downloading LORENE, in the directory Lorene/Codes/Nrotstar. It is a descendant of 
rotstar, a previous LORENE code described in [49] and used in many astrophysical studies (see 
e.g. [56], [108] and references therein). 

nrotstar solves the Einstein equations (3.14)-(3.17) via the self-consistent-field method ex- 
posed in § 3.5.1. For this purpose it uses a multi-domain spectral method developed in [21] and 
[52] (see [53] for a review about spectral methods). At least three domains are employed (cf. 
Fig. B.l): 

• a nucleus, covering the interior of the star; 

• a intermediate shell, covering the strong field region outside the star; 

• an external compactified domain, extending to r = +oo. 

But more intermediate shells can be used, either inside the star or outside it. For instance, the 
nucleus can be dedicated to the fluid interior and a the first shell to the solid crust [110], or two 
domains inside the star can be used to describe a phase transition [111]. 



Not finished yet ! (to he continued...) 
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Figure B.l: Domains used in tlie nrotstar code. 
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